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PREFACE 
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Journal of Aircraft . May - June 1969, pp. 209-214, are copy- 
righted by AIAA, and are used here with their permission. 
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and Mechanical Sciences, Princeton University, in partial ful 
fillment of the requirements for the Ph.D. degree. 
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WIND TUNNEL INTERFERENCE FACTORS 


FOR HIGH-LIFT WINGS IN CLOSED WIND TUNNELS 


Robert Glenn Joppa 


SUMMARY 

A problem associated with the wind tunnel testing of very 
slow flying aircraft is the correction of observed pitching 
moments to free air conditions. The most significant effects 
of such corrections are to be found at moderate downwash an- 
gles typical of the landing approach. 

The wind tunnel walls induce interference velocities at 
the tail different from those induced at the wing, and these 
induced velocities also alter the trajectory of the trailing 
vortex system. The relocated vortex system induces different 
velocities at the tail from those experienced in free air. 

The effect of the relocated vortex and the walls is to cause 
important changes in the measured pitching moments in the wind 
tunnel . 

A method of calculating the interference velocities is 
presented in which the effects of the altered wake location 
is included. The flow fields of a lifting system are calcu- 
lated in free air and in the tunnel, and when compared the 
differences are charged to tunnel wall interference. Itera- 
tive methods are used which require a large computer. The 
tunnel walls are represented by a vortex lattice and the 
results compared with classical methods for the undeflected, 
wake case. 

Results are presented comparing the tail interference 
angles, with and without the effect of vortex wake relocation, 
which show the importance of the wake shift. In some cases 
the tail angle corrections are reduced to zero and may even 
change sign. It is concluded that to correctly calculate the 
interference velocities affecting pitching moments, the 
effects of vortex wake relocation must be included. 



SYMBOLS 


m 

[a] 

Cb} 

b 

b 

g 

c 

c. 


e 

h 


( ) 


n ( ) ( ) 

H 

i/J,k 

L,G 

n 

P 

R ( ) 


( ) 


( ) 


w 


( ) 


Aspect ratio 

Matrix of coefficients of wall vortex elements 

Column matrix of coefficients of wing vortex system 

Wing vortex span 

Wing geometric span 

Wind tunnel cross-section area 

Wing lift coefficient 

Distance downstream to wake roll-up 

Normal distance to a point p from a line contain- 
ing a vortex segment identified by subscript 

Normal distance to a point p from a plane contain- 
ing vortex segments identified by subscript 

Height of wind tunnel 

Unit vectors in the directions X, Y, Z 
Dimensions of rectangular vortex ring (Fig. 8) 

Unit vector normal to vortex ring 
Point having coordinates X, Y, Z 

Vector from point (X,Y,Z) to end of a vortex vector 
S indicated by subscript 

Magnitude of component of vector R, . indicated by 
second subscript ' ' 

Wing area 

Vector representing a vortex segment of strength T 
and length S 

Component of S indicated by subscript 

Unit vector in the direction of the total velocity 
vector at a point 
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Velocity induced at a point 

Velocity component in direction indicated by 
subscript 

Vertical component of wall-induced interference 
velocity 

Width of wind tunnel 

Vector representing a wing bound vortex of strength 

r 

w 

Cartesian coordinate of a point 

Angles defining direction to a point from the end 
of a vortex segment (Fig. 7) 

Circulation strength of a vortex 

Difference between angle of attack in free air 
and in wind tunnel 

Wind tunnel interference factor 

6 Evaluated at tail location 

6 Evaluated at wing location 



I . INTRODUCTION 


The problem of how to do meaningful testing of high lift 
systems in wind tunnels has been with us for some time. That 
wind tunnel testing is necessary for new types of slow flying 
vehicles is evident because the nature of the problems of sta- 
bility and control are different than in flight at cruising 
speeds . 

To obtain the necessary lift at low speed requires that 
incoming air be deflected through a large angle and/or accel- 
erated to a high discharge velocity at a moderate deflection 
angle. In either case the change in angle or increase of 
velocity is no longer small, and so linearized assumptions 
are no longer valid. Pitching moments felt by the airframe 
due to the large turning angle are generally large and non- 
linear, and vary with forward speed as well as with angle of 
httack. 

The gross effects may be estimated by recourse to momen- 
tum methods. Unfortunately, the gross effects are modified 
by real fluid effects that are configuration dependent. Lift 
is developed by real devices such as rotors, fans, and wings 
with flaps. These devices are operated at or near their max- 
imum capability, i.e., near the point of flow separation. In 
many cases, flow separation and re-attachment occur cyclically 
during normal operations, so that linear relationships such as 
between forces and angles of attack, do not usually exist. 

As a result of all this, classical aerodynamic theory, 
which is linearized and limited to small angles, is incapable 
of predicting performance. The only recourse left to the 
designer, then, is to go to the wind tunnel to determine exper- 
imentally the characteristics of a new machine. 

Unfortunately, the wind tunnel introduces its own set of 
problems. While it does indeed permit the solution of the 
detailed problems of separation and mutual interference by 
direct analogy, the quality of that solution depends upon the 
quality of the match of the necessary similarity conditions. 
These are the exactness of the model and the matching of 
Reynolds numbers and Mach numbers. 

High lift systems usually involve rather intricately de- 
tailed parts such as blowing or suction slots, rotors with 
dampered hinges and important elastic properties, or internal 
ducting and fans. The accuracy with which these details can 
be matched imposes some limit on the smallest feasible model 
size; and, in addition, these elements may be the ones most 
sensitive to mismatching of Reynolds number and Mach number. 
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Matching of Reynolds number and Mach number, of course, 
are mutually exclusive except in the case of a full scale 
model. Since the flight speeds of concern are usually low, 
one's first thought is that the test Mach number might be 
increased in favor of a larger Reynolds number, but this is 
not usually possible. At high lift coefficients, local flow 
velocities are often very high and large enough to be affected 
by the local Mach number. Where rotating parts are in use, 
the Mach number of an advancing blade is frequently the con- 
trolling factor. Thus, the test engineer is forced to do 
what he has always done; to accept a lower Reynolds number 
and attempt to extrapolate to full scale results on the basis 
of previous experience. This experience is not extensive at 
present and so he does this very reluctantly, insisting on 
the largest possible model for a given tunnel. 

The wind tunnel also introduces another set of problems 
which are a direct result of the physical presence of the 
boundaries of the test section. The flow from a high lift 

_ sy'stem~has~a— : large-local— down wash— angle-and—' velocity ,__and_in 

free air may require several times its own characteristic 
length to reach final values which may still be very large. 

The wind tunnel walls force the final value of downwash angle 
to be zero and alters both the direction and curvature of the 
flow in the immediate vicinity of the model by an amount which 
is significant with respect to the camber of the lifting sys- 
tem, especially when the model is long (e.g., a rotor, or a 
horizontal tail aft of a wing) . 

That such flow interference exists has of course been 
recognized from the earliest use of wind tunnels, and class- 
ical theory exists for the prediction of the interference 
effects and for the correction of data. Unfortunately, the 
classical work depends on the assumption that the downwash 
velocities are small and that the wake of the lifting system 
goes straight downstream. 

Three methods of coping with this lack of an adequate 
interference prediction theory are available. One can use a 
very small model in available tunnels, build bigger wind tun- 
nels, or develop new theory. A criterion for smallness of 
models was put forth in 1956 (Ref. 1) which suggested that 
the change in curvature of the flow would be sufficiently 
small if the interference angle at the lifting system, cal- 
culated by linear theory, was never larger than 2° . That 
this leads to extremely small models is demonstrated by Fig. 

(1) where it is applied to a helicopter rotor. These small 
models, of course, aggravate an already serious Reynolds num- 
ber problem; and so the industry, still having no adequate 
theory, began in the early 1960's to build larger wind tunnels 
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having test sections of the order of 400 to 1000 square feet 
Even this new generation of wind tunnels is inadequate for 
matching Reynolds number, although the new facilities do per 
mit construction of models large enough that detail can be 
matched with available fabrication techniques. A consider- 
able amount of effort has been devoted to the wall interfer- 
ence problem but a complete solution is still not available. 
This paper is devoted to the development of a new method of 
predicting wind tunnel wall interference for an important 
class of slow flying vehicles. 
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II. DEVELOPMENT AND CURRENT STATE OF 
WALL INTERFERENCE THEORY 

In the classical wind tunnel interference problem, it is 
assumed that the model lifting system can be represented by a 
lifting line and a pair of vortex filaments which trail down- 
stream in a straight, level line from a point near the wing 
tips. A cross-section normal to the flow is examined down- 
stream from the plane of the lifting line, and a pattern of 
other vortex filaments is chosen outside the tunnel walls in 
such a way that the tunnel walls become streamlines of the 
flow. The effect at the model of the added vortices then con- 
stitutes the interference effect of the walls. 

Prandtl presented a solution for the circular wind tunnel 
(Ref. 2) which required only a single pair of vortices outside 
the tunnel wall to cancel, at the wall, the effect of the 
trailing pair inside, but he did not include the effect of the 
lifting line itself. Consequently, his solution is vali d onl y 
at - the plane o f the ~1 Yf ti n g line and cannot give the longitu- 
dinal variation of the interference angles. 

Glauert followed (Ref. 3) with a solution for a rectangu- 
lar tunnel. Since the walls were planes, it was required only 
that each wall become a plane of symmetry of the vortex lines 
inside the tunnel and those outside it, thus leading to a 
doubly infinite set of vortex lines. In the rectangular tun- 
nel there is no problem of how to handle the bound vortex, for 
its external image clearly joins the images of each trailing 
pair. His solution then is valid for points fore and aft of 
the lifting line, and it was possible to show that the effect 
of the tunnel walls was different at the tail than at the wing. 

Other authors have developed solutions for other tunnel 
shapes, but no proper image system has been presented for any 
other shape than the rectangular tunnel. Lotz (Ref. 4) was 
successful in developing solutions for circular and elliptical 
cross section tunnels which accounted for the effect of the 
bound vortex. She added to the image system of Prandtl, a 
potential function expressed in infinite series form, which 
was required to cancel at the wall the normal velocities at 
the wall caused by the bound vortex and also expressed in 
infinite series form. The accuracy of the results depends on 
the evaluation of the truncated series, and no indication is 
given in the original report of the probable error. 

Clearly the basic assumption of the straight downstream 
wake trajectory had to be modified for the consideration of 
the high downwash systems of interest Tiere. The most success- 
ful change to date was made by Heyson (Ref. 5} who let the wake 
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be straight, but at an angle downward until it struck the tun- 
nel floor. The zero size lifting system was represented by a 
point doublet and the wake by a string of such doublets. When 
extending to a finite span wing, a series of such point sys- 
tems are placed side by side; and, since internal singularities 
cancel each other, the result is equivalent to a lifting line 
and a single trailing pair of vortex filaments. The angle of 
descent of the trailing system was taken originally as 1/2 the 
final downwash angle calculated by momentum theory for the 
span-circle mass of air required to produce the lift of the 
system. In a later publication (Ref. 6), he modified this to 
1/4 of the final downwash angle, agreeing with a calculation 
by the author that vortex filaments of a wake move downward 
at approximately 1/5 the final momentum downwash value. Thus, 
the angle of descent used in later work is representative of 
the final wake trajectory, in free air, of the trailing vortex 
system. Image systems are then constructed outside the tunnel 
(rectangular cross-section) . At the point where the trailing 
wake strikes the floor, it is met by the first image wake, 
and they are assumed to change direction and move aft together 
in the plane of the floor. 

With the image system constructed as described, it was 
possible to sum the interference velocities at the model due 
to the external vortex system. It should be noted that the 
doublets, normal to the plane of the downward trailing pair, 
have fore and aft components as well as vertical components; 
and, consequently, longitudinal as well as vertical interfer- 
ence velocities exist. At the floor intersection, only the 
vertical components are canceled; the longitudinal components 
add and are retained. 

Some controversy exists about the degree to which these 
interference calculations are applicable. Evidence has been 
presented (Ref. 6,7,8) to show that good results are achieved 
when calculating interference velocities at the model and 
using them to correct lift and drag. The method has not been 
uniformly successful in correcting pitching moments, however. 

As an indication of the controversy, it may be said that 
another laboratory has offered evidence that wind tunnel and 
flight stability data may agree more closely when no correc- 
tions whatever are applied (Ref. 9) . 

The solutions of Heyson, and others who have tried to do 
something . similar , are deficient in at least two respects. 

The first and most obvious is that the assumed wake position 
is not correct. Others have attempted to improve on the wake 
trajectory by using other assumptions or by modeling experi- 
mentally measured wakes, and then using Heyson ' s computations 
to calculate the interference velocities due to images of 
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these more correct wakes. Results are reported to be little 
changed at the model location, but they are still inadequate 
for pitching moments. 

The second deficiency is the one which is the more impor- 
tant and which no one has yet attempted to account for. This 
is the direct effect on the model of the fact that the wake 
trails along a different trajectory in the tunnel than in 
free air. The effect arises this way. The presence of the 
boundaries (as made evident by the image system) causes upwash 
velocities which are felt everywhere in the tunnel; by the 
model tail and also by the vortex wake itself. The result of 
these upwash velocities is to cause the vortex wake to be 
higher in the tunnel than in free air. This new higher posi- 
tion is different with respect to the tail. For example, if 
the tail is above the wake in free air, the wake will now be 
raised closer to the tail and will induce on the tail a strong- 
er downwash than in free air. This effect may equal or exceed 
the wall or image induced upwash, and thereby dominate the 
pitching moment interference. 
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III. A NEW APPROACH TO INTERFERENCE CALCULATIONS 

A new approach to the problem is offered in this paper 
which attempts to remove the two deficiencies of former meth- 
ods. The interference must be computed for the correct wake 
shape, and the direct effects of the relocated wake must be 
included. In order to do this, the flow field of the lifting 
system must be predicted both in the free air case and in the 
wind tunnel, and the differences in flow velocities be charged 
to wall interference. In order to develop the method, certain 
restrictions to the problem were defined for practical reasons. 

The principal effect which it is desired to show is that 
the relocation of the wake by the interference of the walls 
contributes a major influence on pitching moment interference, 
which may be added to or subtracted from the usual interfer- 
ence calculations. It is not difficult to show that the effect 
of a shift in the wake position will have a maximum effect when 
the wake is only moderately deflected with respect to the tail 
or the plane of a rotor. Figure (2) shows a section taken 
(Trefftz plane) at a location representative of a tail with a 
pair of trailing vortices at a distance h below the tail. 

The downwash is given by the Biot-Savart equation, and is 


w 


TT | [1 + 


b/2 


) 2 ] 


The ratio of the downwash velocity to that experienced when 
the wake is at the same height as the tail, (h = 0), is given 
by 


w 

W (h = 0) 


1 


1 + ( 


h 

b/2 


2 


The maximum rate of change of downwash with height occurs when 



If the length of the model is of the same order as the 
span, and the model is in a level attitude, then this corres- 
ponds roughly to a downwash angle of the vortex wake of about 
16°. Helmbold (Ref. 10), has shown that the maximum lift 
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possible due to circulation alone will produce a wake trajec- 
tory angle of just over 20°. Therefore, the attainable values 
of circulation lift place the wake in the region where changes 
in its location will produce the maximum effect on the down- 
wash at the tail. 

Greater wake trajectory angles are of course produced by 
highly powered lifting systems where the power is used to 
increase the mass rate of flow through the system. Analysis 
of highly powered systems is not included here for two prin- 
cipal reasons. First, the larger downwash angles remove the 
wake vorticity further from the tail plane, and so the effects 
of wake relocation become less important. If the downwash 
angles are large enough, the tail is almost unaffected by 
changes in wake location, and in this case the methods of 
Heyson become appropriate, and indeed have given good results. 

A more practical reason for avoiding larger downwash 
angles is that at some point interaction with the tunnel walls 
produces an impossible situation. In the limiting case of 
hovering inside a test section, the forces measured are clear- 
ly different from those in free air because of recirculation 
of the air. For a range of forward speeds above hovering, 
recirculation still exists in the tunnel where it will not in 
free flight, even near the ground. At speeds just above re- 
circulation, experiments by Rae (Ref. 11) indicate that forces 
measured are so far from what is expected that test results 
are highly doubtful and may be useless. Apparently the rotor 
wash is interacting with the entire tunnel flow and producing 
a large circulation very close downstream in a way which has 
yet to be satisfactorily explained. His test results show 
that a fairly definite point can be determined at which this 
effect (which he calls flow breakdown) disappears and one 
expects credible results. This limit probably determines the 
lower speed bound (maximum downwash angle) for corrections of 
any type. Consequently, this region will not be examined here, 
and the problem will be confined to lifting systems which can 
be said to produce only circulation lift. 

This type of system is simply represented as a lifting 
vortex line with a single trailing pair of vortices. Such a 
mathematical model could represent a simple wing with some 
sort of boundary layer control so that the large values of 
circulation can be developed. It may also represent a heli- 
copter rotor operating in the translational lift region. Since 
we are primarily concerned with the flow field at a distance 
from the model (at the tunnel walls), details near the model 
are of lesser interest and a relatively simple model represen- 
tation can be used. 
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It is assumed that the trailing sheet of vorticity rolls 
up immediately into a cylindrical core of vorticity which can 
be represented by a single filament located at the center of 
gravity of the original vortex sheet. Actually, this assump- 
tion is not really necessary. It only need be shown that the 
effect of the singular representation of one half of the trail- 
ing sheet on the center of gravity of the other half is not 
significantly different from the effect of the real sheet. 

It is demonstrated in Appendix A that the effect of the unde- 
flected sheet trailing from one half of an elliptically loaded 
wing is only 2%% larger than the corresponding effect of a sin- 
gularity at the center of gravity. After roll-up, the vortex 
sheet becomes axially symmetrical and it is easily shown that 
the effect at any external point of a uniform cylindrical vor- 
tex sheet is identical to that of a filament at its center 
having the same total strength. 

Evidence that the wake does roll up quickly is given by 
Sprieter and Sacks (Ref. 12) who report the roll-up distance 
as a fraction of the geometric wing span to be 



0.28 ( 


7R 


) 


In the high-lift case of interest here, fll/C L is about 1.0, 
so the roll-up distance would be of the order of a chord 
length downstream. 

That a helicopter rotor can be represented by the lifting 
line and trailing pair is graphically shown by data taken by 
Heyson, (Ref. 13). Figure (3), taken from NACA TR 1319. shows 
that for a rotor having a momentum downwash angle of 15 , two 
clearly defined vortex cores are already well developed at a 
plane only just downstream of the rotor trailing edge. It 
also shows that the cores are deflected less than one half as 
much as the air mass, calculated by momentum theory. 

In summary, the problem that will be presented is the cal- 
culation of the interference due to the walls of a closed test 
section wind tunnel, on a high-lift wing having a moderately 
large downwash angle, taking account of the direct effect of 
the relocation of the vortex wake on the longitudinal distri- 
bution of downwash. The problem is approached by first cal- 
culating the trajectory of the wake of a simple lifting sys- 
tem and its flow field in free air. The lifting system is 
then placed in a wind tunnel and its new trajectory and flow 
field are compared at the same values of remote wind speed and 
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model circulation strength; differences are interpreted in 
terms of tunnel wall interference. In order to determine the 
flow field in the wind tunnel, a new method of representing 
the wind tunnel walls was developed and is also presented. 
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IV. THE FREE AIR TRAJECTORY 


Figure (4) shows a sketch of the vortex wake representing 
a plane elliptical wing and indicates the induced velocity due 
to an element of the vortex acting at an arbitrary point. The 
element of induced velocity is evaluated by the Biot-Savart 
law, and when integrated over the entire wake, the direction 
of the flow at a point can be determined. The flow direction 
is first determined along an initially assumed wake trajectory 
and the wake is then deflected to assume the calculated direc- 
tion. With the wake now deflected, a new calculation of flow 
direction is made and the solution converges after several 
iterations . 

To facilitate the solution, the vortex system is broken 
into a series of short straight line segments. The bound vor- 
tex lies on the quarter chord line and has a span of tt/4 
times the geometric span, which is appropriate for represent- 
ing an elliptical wing. The first trailing segments lie in 
the plane of the wing, extending from the bound vortex tips 
to the trailing edge. The downstream vortices are assumed to 
spring from the trailing edge at that point and are divided 
into segments whose length is approximately 1/10 of the vortex 
span. The angle of the first segment, being in the plane of 
the wing, is determined by adding the induced angle of attack 
and the effective angle of attack at the plane of symmetry. 

The induced angle of attack of the wing is computed at the 
lifting line by summing the induced velocities of all the 
trailing segments and adding them vectorially to the remote 
velocity. The effective angle of attack is determined by 
assuming two dimensional flow at the plane of symmetry and 
setting the normal component of the local velocity vector 
equal and opposite to the velocity induced by the bound vor- 
tex at the three-quarter chord point. See Figure (5) . 

The direction of each downstream element, in turn, is 
calculated by summing the individual velocities due to all 
other elements at its own upstream end. This direction is 
used to determine the coordinates of the downstream end of 
the segment; the entire string of segments downstream from 
that point is translated so that it stays attached, and the 
next segment direction is determined. Thus, the wake is 
moved into place by sweeping along its length from the wing 
aft in several iterations. 

When a vortex line lies in a plane and follows a path of 
varying curvature, it induces on itself velocities normal to 
the original plane which vary with the curvature. The fila- 
ment, which leaves the wing at a fixed location, curves upward 
from its angle of departure, and so each downstream section 
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experiences an inward deflection from its own upstream ele- 
ments. This vanishes as the trajectory straightens out, but 
it must leave the final straight wake at a smaller vortex 
span than it had on leaving the wing. The iteration process 
must then allow for this lateral freedom, as well as for the 
vertical motion of the wake. 

When the above described process was first attempted, 
simultaneously calculating both downward and inward deflec- 
tions, the computation became unstable after only a few iter- 
ations. This instability was avoided by a double iteration 
process. First, one pass is made calculating only downward 
deflections, and then a second is made allowing only horizon- 
tal or inward deflections. By this stepwise process, a tra- 
jectory can be found which converges after only three or 
four such double passes, and which converges before instabil- 
ity develops. 

It should be noted that the vortex line is physically 
unstable in that curvature of the line causes more self- 
induced curvature. A pair of vortex lines, if disturbed, will 
break up into segments and eventually produce vortex rings. 

An example may be observed in the contrails of jet aircraft, 
where the engine exhaust is drawn into and makes visible the 
cores of the trailing vortex pair. This instability could be 
accentuated by round-off errors in the computing machine and 
places a limit on the number of times an iteration can be car- 
ried out. 

A computer program with instructions and card listing 
for the solution for the vortex trajectory from a lifting 
wing is given in Appendix B. 
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V. REPRESENTATION OF THE WIND TUNNEL WALLS 

While the image systems described earlier are correct, 
and could be used with proper modification for finding the 
interference velocities due to the tunnel walls, they still 
leave something to be desired. Since the vortex wake of the 
lifting system in the tunnel will be curved, the external 
images would also have to be curved; and furthermore, since 
the final solution will have to be iterative, the geometry of 
the image system will have to change also for each iteration. 
These problems can be handled by a computer, but the method 
has some more basic restrictions. Proper images are avail- 
able only for rectangular tunnels and the concept of an image 
implies that the tunnel is of infinite length. Tunnels in 
use for high lift testing are not all rectangular and, more 
important, many of the special tunnels being built today have 
such short test sections that some doubt exists about their 
adequacy. Therefore, in an effort to satisfy these objections 
a new approach was developed. 

In this method the image concept was abandoned and the 
tunnel walls are represented by a vortex lattice. The strength 
of each element of the lattice is found by simultaneously 
requiring that the normal component of velocity vanish at a 
control point in the center of each lattice element. This 
method has the computational advantage that the geometry of 
this system is unchanged during each iteration, and that the 
large matrix of coefficients need be inverted only once for a 
series of computations. 

Further, it is applicable to any tunnel cross section to 
the extent that it can be approximated by a polygon of equal 
length elements, and the effects of finite length can be 
explored. In order to prove the method, it was first applied 
to the classical problem of the undeflected wake. The devel- 
opment follows. 


Problem Statement 

The problem is to find that distribution of vorticity 
lying in the tunnel walls which will prevent any flow through 
the wall due to the action of a lifting system in the wind 
tunnel. The lifting surface is assumed to be uniformly loaded 
and is represented by a simple horseshoe vortex with the trail- 
ing pair undeflected. In principle, any desired distribution 
of lift could be built up of several such simple elements. 

The walls are represented by a tubular vortex sheet of 
finite length composed of a network of circumferential and 
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longitudinal vortices having equal spacing. Helmholtz' theo- 
rem that a vortex filament can neither end nor begin in the 
flow is satisfied most readily by constructing the network of 
square vortex rings lying wholly within the plane of the walls. 
Each square has a vortex strength Fj_ , and each side is co- 
incident with the side of the neighboring square. Thus, the 
strength of any segment is the difference between the strengths 
of the two adjoining squares. The boundary condition that the 
wall must be impervious to flow is satisfied at a control 
point in the center of each square. This results in a set of 
simultaneous equations, one written for each control point, in 
which the unknowns are the F^ . 

A large number of equations results if the tube is very 
long, thus some judgment is required in choosing the geometric 
arrangement. The use of square vortex rings requires a tunnel 
of constant cross-section. One notes that for a wing mounted 
in the center of the tunnel, lateral symmetry always exists; 
and, if the wake is undeflected, vertical symmetry also exists, 
thus reducing the number of unknowns. The trailing edge of a 
finite length tube which represents the long tunnel requires a 
slightly different treatment. At a far downstream section, 
only longitudinal vorticity should exist. This is represented 
by elongating the last ring of squares by a large amount, while 
keeping the control point at the same location with respect to 
the last circumferential station. Figure (6) shows the arrange 
ment for a rectangular tunnel with filleted corners. 

Equation Setup and Solution 

A right-hand axis system is established with the X-axis on 
the longitudinal centerline of the tunnel, positive downstream. 
The Y-axis is taken positive upward and the Z-axis positive to 
the right side of the tube facing downstream. 

Since the surface of the tunnel is to be made of square 
elements, its cross-section is a polygon of equal segments 
arranged to approximate the desired cross-section shape. In 
this development, the cross-section will be assumed to be 
symmetrical about the X , Y plane. 

In general, the velocity induced at any point p (Fig. 7) 
due to a vortex segment may be written: 

V = (cos + cos p 2 ) v (1) 
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where v is a unit vector to establish direction. The terms 
required are written as follows: 
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Finally, the velocity induced at a point due to a vortex seg- 
ment is: 
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One could then add the contributions of all four sides 
Of a vortex square, but it is more convenient to take advan- 
tage of the lateral symmetry and sum the effects due to a 
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pair of symmetrically located vortex squares of the same 
strength. The arrangement is shown in Fig. (8) and the 
following equation results: 
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Similarly, the velocity induced at point p by a simple 
horseshoe vortex located in the center of the tunnel is de- 
rived from Fig. (9) using Eq. (1). Summing the contributions 
from the three segments yields : 
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iff 


The normal velocity at a point on the wall is constructed 
by taking the dot product of the induced velocity vector with 
the unit outer normal at that point. V n = V • "n . The nor- 
mal is constructed using the cross product of a unit vector 
in the downstream direction and a vortex ring vector lying in 
the Y - Z plane 
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| i x (R 1 -R 2 ) j 


The boundary condition is expressed at each control point by 
summing all the normal velocities due to the wall vortex rings 
and setting it equal and opposite to the normal velocity 
induced at the same point by the wing vortex. The result is 
expressed in a matrix equation 


[A] 
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in which the {T} are the unknown strengths of the wall vor- 
tex elements, and the matrix [A] is fixed by the dimensions 
and shape of the tunnel and the locations of the vortex rings 
and control points. The column {B} describes the influence 
of the lifting wing at the tunnel walls, and is developed from 
the dot product of Eq. (4) with the unit outer normal at each 
control point. 

Because of the lateral symmetry assumed in writing Eq. 

(3), it is necessary only to take control points on one side 
of the tunnel. If the wing is also placed on the vertical 
and the tunnel is vertically symmetrical, then the Ti 
will also be symmetrical but of opposite sign. It is then 
necessary only to take control points in one quarter of the 
tunnel. The matrix [A] is inverted, since it is fixed for 
a given tunnel shape, and the values of may then be 

found for a variety of wing spans by changing only the column 
matrix {B} . 

Once the are known, the induced velocity due to the 

walls can be calculated at any point in the tunnel by the use 
of Eq. (3) summed over all the vortex rings in the tunnel 
walls. The interference is expressed as an angle whose tan- 
gent is the vertical component of interference velocity, w , 
divided by the tunnel wind speed, V . In the linear, unde- 
flected wake case, the tangent is approximately equal to the 
angle. Results are expressed in terms of the classical inter- 
ference factor 6 , defined by the equation: 



The factor is computed in terms of wing circulation and vor- 
tex span 


5 


_w C_ 
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w 


Results are presented graphically to show the longitudinal 
variation of the factor 6 for different wing spans in a 
variety of tunnels. A computer program with instructions and 
card listing for the solution of the interference factor 6 
is given in Appendix C. 
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Results and Comparison with Classical Results 

In order to test the validity of the method, it was com- 
pared with classical solutions where those were available. 
Results of calculations made for three representative tunnel 
shapes are presented in the form of graphs of the wall inter- 
ference factor § . Values of 6 were calculated at points 
along the tunnel centerline from the wing location downstream 
for several values of wing vortex span. These are presented 
for a circular, a square, and a 3:5 rectangular tunnel in 
Figs. (10), (11), and (12). The average value of this inter- 

ference factor over the vortex span of the uniformly loaded 
wing was also calculated and is shown as a function of vor- 
tex span for each of these tunnels along with the centerline 
values in Fig. (13). 

Square tunnel . — Glauert's concept of an infinite array 
of images of the wing located outside the tunnel is applicable 
only to rectangular (including square) tunnels and has been 
applied by Silverstein and White in Ref. (14). Results are 
presented there for square and 2:1 rectangular tunnels; only 
the square tunnel results are used here for comparison, 
since 2:1 tunnels are not common. 

The number of line segments, each corresponding to the 
side of a vortex square, to be used to adequately represent 
the square tunnel cross-section was determined by making a 
series of calculations with increasing numbers of segments. 
Fig. (14) shows the results of using 12, 16, and 20 segments 
to make up the periphery of the square cross-section. The 
results for 16 and 20 segments differ only slightly and 
correspond very closely to the data taken from Ref. (14) . 

The excellent agreement shown indicates that 16 segments are 
enough to represent satisfactorily the square cross-section 
tunnel . 

Circular tunnel . — In the case of the circular tunnel, 
no exact solution is available for the downstream interfer- 
ence factors, so two approximate results are compared with the 
new calculations in Fig. (15) . The results presented by Lotz 
(Ref. 4) depend on the value of a truncated infinite series, 
and the reference gives no indication of the accuracy expec- 
ted in its evaluation. The result taken from Silverstein and 
White (Ref. 14) was arrived at by following their suggestion 
that the downstream interference factors for the circular 
tunnel be taken as the same as for the square tunnel of the 
same area. 

Four different approximations to the circular tunnel were 
used for this calculation. Two regular polygons having 12 or 
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16 sides were used for the cross-section shape; each was rota- 
ted so that either points or flats of the polygon were at the 
top and side centerline. All four calculations yielded the 
same curve, with values within one-tenth of one percent. Thus, 
it is concluded that a twelve-sided polygon is adequate to 
represent the circular tunnel. 

Length effect . — The effect of length of the tunnel to 
be used in calculations was explored for the circular tunnel. 

A twelve-sided polygon was used in the calculation, with the 
model vortex span equal to 0.4 of the tunnel diameter. It is 
evident from Fig. (16) that a length- to-diameter ratio of 3 
or 4 is ample for convergence. The reason for this may be 
seen in an examination of the distribution of the wall vor- 
ticity. The bound vortex of the wing requires some circum- 
ferential vorticity in the walls, but only in the region quite 
near to the wing. Longitudinal vorticity is not required far 
upstream, and far downstream only longitudinal filaments exist 
to control the trailing pair from the wing. By using the 
artifice of a very long last ring, the proper conditions are 
met far downstream, and the vortex lattice need only be long 
enough to provide the circumferential vorticity needed in the 
immediate vicinity of the wing. In fact, all the vorticity 
in the circumferential rings is quickly transferred to the 
longitudinal filaments. 

Figure (17) shows the wall vortex strengths taken from 
calculations made for circular tunnels of various lengths. 

The circumferential vorticity strengths were taken at the 
floor near the center of the tunnel where they are the strong- 
est; the longitudinal vortex filament strength is that along 
the side wall at model height. It is evident that the details 
of the distribution are not strongly affected by the presence 
or absence of tunnel walls more than about one diameter up or 
downstream from the wing. 


Conclusion 

The excellent agreement shown by the examples presented 
verifies the hypothesis that the walls of the tunnel may be 
adequately represented by a rather coarse network of vortex 
rings. The advantage of this method is that any tunnel cross- 
section can be represented by using an equivalent polygon of 
16 or more equal length sides arranged to approximate the 
actual geometry. 


23 



VI. THE FINAL SOLUTION 


The solution for the wake trajectory in the wind tunnel 
is an iterative combination of the free air trajectory solu- 
tion and the wind tunnel wall vortex lattice solution. The 
lifting system, represented by a horseshoe vortex, is placed 
inside a vortex lattice tube representing the tunnel, and is 
given an initial value of circulation strength and an unde- 
flected wake. A solution is found for the wall vorticity 
exactly as described in the earlier section. The wake loca- 
tion is then found exactly as in the free air solution, with 
the exception that the velocities induced by the wall vortic- 
ity found for the undeflected wake are added to those induced 
by the wing on itself. After an equilibrium trajectory is 
found, a second solution for the wall vorticity is made with 
the wake in its deflected position, followed by a second iter- 
ation of the wake location. In general, the two systems do 
not interact strongly for the short span to tunnel size ratios 
one expects to use in testing of high lift systems; and so 
only two or three such cycles are usually necessary for con- 
vergence. 


Determination of the Interference Factors 

In order to find the total interference effect, one 
should compare the flow patterns of the system, operating at 
the same conditions, in and out of the tunnel. The same con- 
ditions, as used here, mean at the same circulation and remote 
velocity. When the solutions are complete, they yiel the 
complete velocity field both in free air and in the c;:;. .:J1, 
as well as the separate contributions to that field by the 
wall vortex lattice and the lifting system. 

The interference velocities are then defined by stating 
that the difference between the velocity at a point in the 
tunnel and the velocity at the same point in free air is the 
total interference velocity. Both the horizontal and vertical 
components of the interference velocity should properly be 
considered, but because the moderate wake deflections of the 
examples considered here cause only very small longitudinal 
interference (3% in the extreme cases) , only the effects of 
the vertical component are presented. The vertical component 
of the interference is felt as a change in the angle of attack 
so it is convenient to present the interference in those terms. 
Thus 


Aa a tunnel “ a free air 
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These angles are not small enough to allow the use of 
the small angle approximation so they are defined by their 
tangents. 



This data is usually presented in terms of a value of 6 
defined by the equation 

- 6 if C l 


but since we are comparing at equal values of T instead of 
, we use the relation 
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A computer program listing is given in Appendix D for 
the combined solution for the interference factor 6 for a 
lifting wing with deflected wake in a closed tunnel. 


Resul ts 

Calculations are presented for a plane wing, at lift 
coefficients approaching the maximum theoretically possible 
for an unpowered system. In order to achieve the highest wake 
deflection angles, the aspect ratio of sample calculations was 
taken at 3.0 so that high C^/IR values could be attained. 
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The wing vortex span was taken as one half the tunnel width, 
and the tunnel had a rectangular test section of height to 
width ratio 1:1.5 . 


Figure (18) shows the trajectory of the wake in free air 
and in the wind tunnel for the sample wing. The differenced in 
location of the wake in the tunnel is evident. In Fig. (19) 
the value of the interference factor 6 is shown as a func- ’ 
tion of C T /Sl at the location of the wing and for three tail 

I-l „ „ -5 


locations assumed to he on the tunnel centerline. 


The tail interference angle is taken as the difference 
between the interference angles at the wing and at the tail, 
and presented as the difference between the values of 6 at 
these two points. Figure (19) also shows the tail interfer- r 
ence factor (6 - 6 ) . This curve shows that, for the geom- 

U W 

etry chosen, the pitching moment corrections may become small 
or even negative at the higher lift coefficients. 

In order to demonstrate the effect of the wake shift, ; 
Fig. (20) was prepared for comparison with Fig. (19)'. The 
same factors were calculated, but the contribution of the 
deflected wake was left out. The interference angle was' cal- 
culated using only the velocities induced by the wall vortex, 
lattice. The wake location as computed in the tunnel was used, 
so these results accurately represent interference- velocities 
based upon only the wall induced effects. Figure (20) also 
shows the tail interference factors calculated using only the 
wall induced velocities. The importance of including the : 
direct effects of the wake ' relocation is shown when Fig. (20) 
is compared with Fig. (19). • 

Tail location is an important parameter, for if the tail 
is initially below the vortex wake in free air, then the wake 
shift upward in the tunnel will accentuate the wall induced 
upwash. Figures (21) and (22) show this effect for tail 
heights of 0.2 and 0.4 times the vortex span below the wing, 
as well as the reversal which takes place when the wake moves j 
past the tail location. 

In the preceding examples the interference angle factors 
were calculated at fixed locations in the tunnel, and do not 
necessarily represent a physically realizable vehicle. The 
results can be interpreted to represent a tilt-wing type 
vehicle in which the body is constrained to a constant angle 
of attack. 

For the case where body attitude changes, it is necessary 
to calculate and compare flow angles at the tail in freie air' 


26 



with those in thie tunnel at angles of attack appropriate for 
the same wing circulation. An example is presented in Fig. 

(23) for a case where wing and tail are fixed to a body and 
rotate as a unit. The tail is located above the plane of the 
wing (0.2 of the vortex span) and three tail lengths are 
shown. The interference factor shows a minimum where the tail 
passes through the height of the vortex wake. The large varia- 
tions of the factor indicate the importance of accounting for 
the wake shift and for actual tail position. 
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VII. DISCUSSION OF RESULTS 

In this section the results and their implications will 
be discussed in some detail. Some examples will be worked 
out showing how corrections would be made using these inter- 
ference calculations, some of the difficulties encountered in 
making corrections, and how these difficulties may be resolved 
by modifying the test program. Additional discussion considers 
the adequacy of the mathematical model, computational problems, 
and suggestions for possible future modification or growth of 
this method. 


Examples of Corrections of Test Data 

The results presented in the previous section are in the 
form of the factors 6 W used to calculate the correction to 

the angle of attack at the wing, and (6*. - 6 W ) used to cal- 
culate the difference in angle of attack at the tail from that 
at the wing. These values will be used here to compute exam- 
ples of actual corrections that should be applied and show 
their effects on final data. 

The factor 6 W is used to calculate the interference 
angle at the wing in the following formula 

Aa = 5 w C °L 


where Aa is the increase in angle of attack at the wing 
caused by the restriction of downwash by the tunnel boundaries. 
For the examples presented earlier, the following values 
result. The wing has SI = 3 and its vortex span is one-half 
of the tunnel span. The wing area to tunnel cross section area 
ratio is then 2 /tt^ , assuming a vortex span ratio of tt/4 . 
From Fig. (19) , the value of 6 W is almost constant at the 
wing up to C^/Sl =0.5 and is only changed by 10% out to 

C L /iR approaching 1.0. The table shows values of the angle 
of attack interference at selected lift coefficients. 
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deflected wake 

straight wake 

V* 

C L 

6 

Aa 

deg 

iC ° t 

Aa 

deg 


0.0 

0.0 

0.111 

0.0 

0.0 

0.0 

0.0 

0.5 

1.5 

0.115 

2.01 

0.0525 

1.94 

0.0507 

0.7 

2.1 

0.120 

2.93 

0.1076 

2.72 

0.0995 

0.9 

2.7 

0.130 

4.08 

0.1925 

3.48 

0.1642 


The Aa shown is a correction to be added to the angle of 
attack measured in the tunnel. In free air the wing would 
have to be at the higher angle in order to produce the same 
lift as in the tunnel. 

When the angle of attack is corrected the lift vector is 
rotated by the same amount. The effect of the rotation of the 
lift vector then causes a component of the lift to appear as 
an additional drag, the magnitude being equal to the lift co- 
efficient multiplied by the interference angle in radians. 

This result is also shown in the table above. 

If the wake was not deflected, the value of 6 W would 
be constant at all lift coefficients, and the corrections 
would have been smaller. The corresponding values of Aa 
and AC Dt for the undeflected wake are also shown in the 

table. Comparison of the corrections shows that only small 
changes, of the order of 15% of the drag correction, are due 
to wake shift. Since the total drag correction is of the 
order of 25% of the induced drag at the highest lift coef- 
ficient, this change is less than 4% of the measured dra . 

Calculating the difference in interference at the tail 
shows a more dramatic effect. In the normal case (undeflected 
wake and low cV/® ) where 6t and 5 W are constant over the 
range of C L of interest, one calculates the difference in 
angle of attack at the tail and the wing caused by the inter- 
ference and uses this angle to calculate a correction to the 
pitching moment. Since the tail experiences a greater 
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interference angle than the wing, the^ moments measured in the 
tunnel are more negative for positive lift coefficients. 
Because the interference angles are proportional to , the 

effect is to measure a larger negative value of the slope 
d C^/d in the tunnel, making the model appear more ' stable 

than it would be in free air. o- : - 


Because of the wake de flection , v -,- the tail; angle ; correction - 
will be different from what it .would, .be swi thou t' wake deflec- 
tion. The curves of Fig. (19), (21), and (23) show this for 

three dif ferent^ examples-;. :,- > ciar.zu v:*.- :0. cz - : v •: 7. . 

To calculate the change in pitching ^moment requires know- ; 
ledge of the characteristics of the horizontal tail. For an 
example calculation ; let :us - assume; /that the ; tail length t is equal 
to the vortex span, the tail volume coefficient = 1.0 , the 

tail- aspect ratio is about the same as n tffe' wing^' and lias a lift 
curve slope of n/radian. Then the correction to the pitching 
moment would be -.zvz .c : .ns -..tr. >-• 
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The following table compares the corrections for the 
several cases with those expected when the wake goes straight 
back and the tail is at wing height. In the tilt wing case, 
the tail remains fixed at wing height while the wing rotates 
to increase lift. The column headed low tail'is also a tilt 
wing, but the tail is fixed in the tunnel at 0.2b below the 
wing height. In the moving tail case, the tail is assumed 
attached to the wing at 0.2b above- the plane of the wing, and 
moves as the wing rotates in the tunnel. 


C L 

straight wake 
ac m 

tilt wing 

ac m 

low tail 
ac m 

moving tail 
ac m 

0.0 

0.0 

0. 0 . . .. 

0.0 

o 

• 

o 

0.9 

0. 0636 - . . 

0.0522. 

0.805 

0.062 

1.5 

0.105 

0.0679 

0.1482 

0.134 

2.1 

0.147 

- 0.0535 o 

. 0.209, 

0.268 

2.7 

0.189; .,-.-: 

0.0 . 

0.2325 ... 



The tabulated values are plotted in Fig. (24)' to show the 
correction to the pitching moment coefficient for the sev- 
eral cases. If the wake is not deflected, the interference 
would be proportional to C L as shown, and the apparent 
interference is just a change in the stability derivative, 
d C M / d C T , of the aircraft. For the case shown this amounts 


to a change in that derivative of 


dC 


M 


dC. 


= 0.07 and is inter- 


preted as a change in the location of the center of gravity 
for neutral stability of 1 % of the wing mean aerodynamic chord. 


The other cases are not as simple. The effect of the wake 
shift changes the correction very much and how it does so is a 
function of the exact location of the tail with respect to the 
wing. For the case where the wing tilts and the tail stays 
fixed in the tunnel at the height of the wing, the total inter- 
ference may be seen to be the same as for the undeflected wake 
at low C L , but reach a maximum and decline to zero at high 
Cj, . If the tail is lower than the wing, the wake shift effect 
causes the interference to be larger than in the undeflected 
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case because the wake moves closer to the tail. In the case 
where the entire aircraft rotates so that the tail starts 
above the wake and moves past it, the curve shows a reversal 
of initial trend and finally deviates very markedly from the 
no-deflection case. 

The tilt-wing case is perhaps the most interesting of the 
three cases. At low C L values, the corrections are identical 
to those for the undeflected wake, and the stability level in 

dC 

the tunnel is apparently too high by A -. 7 7 - = - 0.07 . At 

ac L 

about Cl =1.5 , the interference effect is now constant, 
so the apparent stability is the correct value. However, a 
constant AC M is introduced which corresponds to a change in 
stabilizer angle of about 1.24°. At C L = 2 . 7 no correction 
in stabilizer angle will be required, but the apparent stabil- 

dC 

ity is now less than the correct value by A - 0.13 . The 

c* Ij 

effect of this change in pitching moments is to move the loca- 
tion of the neutral point a distance of 20 % of the wing chord 
over the range of available lift coefficients. This is about 
the same as the usual allowable movement of the center of 
gravity of a normal aircraft. 

These three cases taken together show that the fact that 
the wake does move with respect to the tail causes the pitch- 
ing moment interference to vary widely; in the examples, from 
zero to nearly twice the values calculated in the usual way 
assuming no wake deflection and tail fixed on tunnel center- 
line. Because of this wide variation it is not possible to 
generalize on the results beyond saying that the interference 
is dependent on the configuration of the aircraft and the wind 
tunnel, and must be calculated for each case. Because the 
variations of interference are of the same order as the linear 
interference and may be of either sign, they are certainly too 
large to be ignored. 


Difficulties in Application 

Actual application of these interference calculations is 
not as easy as presented above, particularly with respect to 
the computation of the pitching moment correction. As this 
correction was presented earlier, it was presumed that the tail 
effectiveness was represented by the derivative d c M /da t and 

that this value was a constant. In the normal airplane this 
is often so, but in the case of the STOL aircraft one cannot 
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make that assumption. The specific difficulties are that the 
local flow angles may be so large that the lift curve slope 
dC /da, is in a nonlinear range, and that the dynamic pressure 
a t 11 the c tail may not be anywhere near the free stream value due 
either to being immersed in low energy wakes from wing flaps 
or high energy wakes from propulsion devices. Consequently, 
it is usually advisable to measure separately the tail effec- 
tiveness by making several runs at different stabilizer angle 
settings and computing directly from this data the values of 
dC^/dat over the range of lift coefficients of interst. This 
much is often done in ordinary wind tunnel work and is even 
more important in the testing of STOL aircraft. 

An additional consequence of the wake shift is now appar- 
ent. The energy wakes are shifted in position and so are 
likely to change the dynamic pressure at the tail. While the 
process described above of measuring the tail effectiveness 
derivative will allow correction under the conditions of test 
in the wind tunnel, these are different from free air condi- 
tions. What is desired is that the tail in the wind tunnel 
be placed in the same air conditions that it would experience 
in free flight. Since the wake in the tunnel is in a differ- 
ent place than in free air, the tail should be moved to occupy 
the same position with respect to the wake. 

The present method allows one to calculate in advance of 
the test program what the wake shift will be for each value of 
the wing circulation. A model could be constructed so that 
the tail height would be adjustable. Stability testing would 
then be done at several positions of the tail to produce a 
family of curves of pitching moment, each one of which will 
be valid for a given lift coefficient, and final data will be 
a composite curve taking data from the several curves at the 
appropriate points. If the wake shifting of the air impinging 
on the tail is the same as that of the vortex cores, and the 
tail is moved that amount, then the wake shift effect on the 
tail moment correction is reduced to zero and only the wall- 
induced effects would be necessary. Variations of induced 
velocity across the span of a model are not large (of the order 
of 10% or less) for models less than two- thirds of the tunnel 
width, and so this method appears to have promise. 

Another uncertainty in the application of these interfer- 
ence results stems from the estimate of the vortex span and 
the resulting value of the circulation strength which is cal- 
culated using the Kutta-Joukowski law. It is apparent that 
this value should be estimated rather carefully before apply- 
ing interference corrections to the data. It may be desirable 
to make some attempt to measure it directly by locating the 
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vortex trajectory in the tunnel. It should be mentioned in 
passing that this is not a new problem and it has always been 
necessary in applying classical corrections to make this esti- 
mate: because the corrections are larger at higher lift co- 

efficients, the estimate is more important. 

■ .? .• A 

Discussion of - Accuracy and Computation Method 

It will have become apparent in the above discussion that 
the quality of the interference calculation depends on the 
representation of the lifting system and the resulting accu- 
racy of the free air flow fields. It is recognized that, if 
one could actually predict the real flow fields with a high 
degree of accuracy, the wind tunnel would no longer be neces- 
sary; and that, if the accuracy is poor, the interference cal- 
culation will have little value. This statement is not as 
contradictory as it may seem, because there is a difference 
between the detailed effects felt in the near field and the 
gross effects in the far field. Regardless of: how it may be 
produced, lift is a result of the generation of circulation 
about some location fixed in the flow field. Consequently, 
if lift is measured and' the! vortex span carefully estimated 
or measured, the induced effects at points as far away as the 
tunnel walls are very well predicted by the Biot-Savart law. 

A wind, tunnel program is designed, to measure more detailed 
effects, particularly those\ due to local flow separation and 
those due to mutual interference of the components of the air- 
craft on each other. No one at this time realistically expects 
to be able to predict these complex events and so replace the 
wind tunnel with a computer. Since the interference calcula- 
tions presented here depend only on the gioss induced effects, 
the accuracy should be adequate for the purose. The represen- 
tation of the model may be improved as much as desired by super- 
position of additional vortex systems, and should be modified 
for other configurations, but the effects at the tunnel wall, 
and therefore the wall vorticity and the resulting induced' 
velocities, will not be changed very much. What such improve- 
ment and modification wili do is account more accurately for 
the direct effect on pitching moments due to wake shift. Cer- 
tainly such work should be done, but the wide variety of 
arrangements possible preclude any generalization in advance' 
and so it will be done on an ad hoc basis. 

Some remarks are in order on the convergence of the numer- 
ical solution, and the instabilities expected in it. Any dif- 
ficulties to be found would be expected in situations where 
the wake was forced to curve most sharply, and this would be 
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when the wing is inside a tunnel and operating at the highest 
lift coefficients. A detailed study was made of such a tra- 
jectory over seven iterations 1 for the aspect-ratio' 3- wing at 
C l //R about 1 . "Two- regions of the wake were selected which 
exhibited the two areas of concern--- instability and conver- 
gence . - ' '• r - 7 - - ar •" 5- Vir-. ’e 

• •• • *. - 

It was expected that in regions of sharp curvature the 
self-induced effects of adjacent segments of the vortex, made 
somewhat unreal by being broken up into short straight sections 
and aggravated by round-off errors, would initiate local cur- 
vature anomolies and cause the solution to degenerate. 

This effect did indeed appear as a wavy motion of the 
segments alternating around a mean- line. -Two or- three such 
zig-zags appeared in the second and third iterations and about 
twelve segments were involved in the seventh. The amplitude 
of these motions grew slowly and did not reach 20% of the - r 
length of the segments until the seventh iteration. This cor- 
• • re-sponded* : to-'a’ deviation of -the- segment direction of 13° or 
less from a'raean line drawn ‘through* -‘them .' •' ^These- waves' disap- 
peared, in ‘the seventh iteration, at about"' one wingspan down- 
stream from the wirtg where the slope of the trajectory had 
become nearly' constant-. The effects" of* the%e small- changes 
of direction - were judged to be negligible and so no r smoo thing 
- sub-routines were used. - **" 4 

- Convergence was examined at' a point- one vortex span down- 
stream from the wing where the trajectory of the vortex line 
was straight over a length of about one span". The locus of 
points of intersection of the vortex line and the tunnel cross 
section was found to be a spiral over' ! the seven iterations. 
Convergence was approximately logarithmic with each motion 
from one iteration to the next being one-half to one-third of 
the previous one. Thus, the convergence is so rapid that the 
fifth iteration moves the wake" less than 1% of the wingspan. 

One concludes from the above" that the solution is quite 
well behaved and no conflict exists-' between convergence and 
stability. Acceptable convergence" is : had at the fourth iter- 
ation, and the growing instability is still acceptable at- the 
seventh, leaving a wide region of choice for the user. 

Future work could well be done ; on r approximate methods of 
predicting wake deflection; for example by choosing a general 
form for the trajectory' curve, and" finding'" its amplitude at 
only a few points. Certainly other approximations will -sug- 
gest themselves. - - - - ■ ~ " ' 



The results presented here, and the method of approach, 
appear to provide as near to an exact solution as is likely 
to be found, and may be used as a standard to which approxi- 
mate and more convenient methods may be compared. 
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VIII. CONCLUS IONS 


The problem of determining the wind tunnel wall interfer- 
ences for high lift wings or lifting systems for slow flight 
has been examined, and a new method of calculating the inter- 
ference effects has been developed. It has been shown that 
the most significant interference is on the measured pitching 
moments and the apparent longitudinal stability of an aircraft 
having a tail, or at least having a longitudinal characteristic 
dimension of the order of its spanwise dimension. The inter- 
ference is a maximum when the system is operating at moderate 
downwash angles which are attainable with lifting systems using 
only small amounts of power and which can be represented by 
passive systems in potential flow. 

The solution developed is based on the use of a vortex 
lattice to represent the tunnel boundaries, and takes into 
account the direct effect of the interference-caused reloca- 
tion of the vortex wake on the flow direction in the region 
of the tail. A method of testing is proposed which can mini- 
mize this effect. 

The following conclusions may be stated. 

1. Representation of the wind tunnel boundaries by a 
vortex lattice system may be used to calculate 
interference velocities for a tunnel of arbitrary . 
cross-section. 

2. Simplified representations of lifting systems may 
be used. The vortex span and point of origin of 
the trailing system are the most important choices. 

3. Wall induced velocities cause the vortex wake and 
high or low energy wakes to be deflected less in 
the wind tunnel than in free air. 

4. • The relocated vortex and energy wakes cause dif- 

ferent flow angles and velocities to be felt at 
the region of a tail and these effects are prop- 
erly charged to tunnel boundary interference along 
with the wall-induced velocities. 

5. The direct effect of the vortex wake shift on a 
tail may be of the same order as the usual wall- 
induced velocities and may be of either sign. 

6. The amount and direction of wake shift effects 
depends strongly on the tail location and so 
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effects must be calculated for each configura- 
tion of interest. 

7. Wake shift effects may be reduced or avoided by 
testing with models whose tail heights can' be 
adjusted to match the energy and vortex wake 
locations for particular regions of interest. 

8. The numerical calculation presented converges 
rapidly (in about three to four iterations), 
but may develop instabilities if carried beyond 
seven or eight such iterations. 

9. The quality of the solution presented is as near 
an exact solution as practical representation of 
a lifting system will permit, and should serve 
to guide the formation of approximations and as 
a standard to evaluate them. 
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Fig. I Tunnel size required to limit wall interference to 2°. 
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Fig. 6 Representation of a rectangular tunnel with corner fillets by a vortex lattice 
of square vortex rings lying in the tunnel walls. 
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Fig. 9 Definition of distances for a horseshoe vortex representing 
a wing located with its midspan at the origin of coordinates. 
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Fig. 13 Effect of wing span on average interference factor 
and the centerline interference factor at the wing. 
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Fig. 15 Comparison of interference factors with classical 
values for a circular tunnel. 
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TUNNEL LENGTH TO DIAMETER RATIO , L/W 
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Fig. 16 Effect of tunnel length on interference factors for 
a circular tunnel. 
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Fig. 17 Effect of tunnel length on well vorticity distribution 
for a circular tunnel . 
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TAIL INTERFERENCE FACTOR (8f-8 w ) INTERFERENCE FACTOR 8 




Fig. 19 Interference factors at wing and tail including 
wake relocation effects. Tail on tunnel centerline. 
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TAIL INTERFERENCE FACTOR ( S t -S w ) INTERFERENCE FACTOR 8 

INDUCED BY WALLS ONLY INDUCED BY WALLS ONLY 




Fig. 20 Interference factors at wing and tail using only 
wall-induced effects. Tail on tunnel centerline. 
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Fig. 21 Interference factors at wing and tail at .2 b v 


below tunnel centerline. 


61 




62 


Fig. 22 Interference factors at wing and 
tail at A b v below tunnel centerline. 





DOWNWASH FACTOR C L / ft 

Fig. 23 Interference factors at wing and tail. 
Effect of tail displacement included. Tail height 
.2 b v above wing plane. 
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PITCHING MOMENT CORRECTION , AC 



WING LIFT COEFFICIENT , C L 


Fig. 24. Pitching moment corrections for several tail locations. 
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APPENDIX A 


COMPARISON OF THE INDUCED VELOCITY OF A DISTRIBUTED 
VORTEX SHEET WITH THAT DUE TO A SINGULAR VORTEX 


Betz* has shown that the first moment (center of gravity 
location) of a group of vortex filaments in a trailing vor- 
tex sheet is constant as they move about in the process of 
rolling up into a cylindrical arrangement. It is well known 
that the spanwise location of the center of gravity of the 
vortex sheet trailing from an elliptical wing is at tt/4 
times the semispan, measured from the plane of symmetry of 
the wing. It is also well known that the induced velocity 
at some large distance from the vortex sheet may be computed 
accurately by replacing the vortex sheet with a single vortex 
of the same total strength located at the center of gravity 
of the sheet it replaces. What is not widely known is the 
variation close to the sheet when this substitution is made. 
The following analysis is presented to show the ratio of the 
induced velocity in the near field computed using the trail- 
ing sheet, to that computed using a concentrated vortex loca- 
ted at the center of gravity of the sheet. 

Consider the Trefftz plane, but just behind an ellip- 
tically loaded wing, as shown below. 



*Betz, A., "Behavior of Vortex Systems," NACA T.M. 713, 
June 1933. 
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The circulation on the wing is given by 



and the strength of the vortex trailing from the point y is 

ar =(fb 


This element of the vortex sheet induces a downwash velocity 
at a point y 0 



dr 

4TT(y 0 - y) 


These equations are combined, and non-dimensionalized by 

letting y = and Yo = • The integral is evaluated 

only over 0 < y < 1 because we are only interested in the 
effect of one half of the wing on the other half. 


w 


- r c 




j 

i. 


V dy 


4TT f V ° <Yo - y)7w 


7 


The integral can be put into a standard form by making the 
trans forma tion 


Then, 


x = Yo - Y 

Y = Yo ~ x 

2 2 2 
Y = Yo ~ 2 y 0 x + x 

dy = - dx 
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and the limits of integration become 

when y = 0 , x = y 0 
when. y = l ,x = y 0 -l 

Then 

w = - r 2_ J y ° " 1 (y° ~ x) dx 

y °y 4 tt j y 0 Xy/ (1 - yo 2 ) + 2y 0 x - x 2 


This is integrated for values of - 1 < y 0 < 0 , using 
integrals number 161 and 182 from Pierce, A Short Table of 
Integrals . Ginn and Company, 1929. The result is 


To 

IT Y9 . - In 

- yo \] 

y °y 4rr ^ 

L 2 /rr^ ln . 

i-y» 2 J 


Now compare this solution with that of the simpler case, where 
the total circulation, - T 0 , is assumed to be concentrated 
at y 0 = rr/4 • b/2 , and find its effect on the other side of 
the wing. We have, then 


-IV 


w 


yo, 


4tt 


b / _yo__ n\ 

2 v b/2 4 / 


The ratio of the downwash due 
single vortex is 


to the sheet to that due to the 


R = 


w (sheet) 
yo,. 


w 


y 0y (single) 


/ i__X_ 

V 4 b/2 


)[*- 


yo 


In 




kr^)] 


We are particularly interested in the value when y 0 = tt/4 , 
and that value is 

R = 1.02566 
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The graph following shows the variation of this ratio over 
a range of distances from the wing. 



X 

b/2 



X 

b/2 


DOWNWASH ALONG THE EXTENDED LIFTING LINE 

R is the ratio of downwash due to a vortex sheet 
trailing from one half of an elliptically loaded wing to 
the downwash due to a single trailing vortex of the same 
strength located at the center of gravity of the trailing 
sheet. 
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APPENDIX 6 

PROGRAM TO COMPUTE THE WAKE TRAJECTORY 
OF A VORTEX PAIR TRAILING FROM A FINITE WING 


PROGRAM FRAIR (INPUT, OUTPUT, PUNCH ,TAPE5=INPUT,TAPE6=0UTPUT, B 

1TAPE7=PUNCH) 9 

8 

THIS PROGRAM IS WRITTEN IN FORTRAN IV FOR THE CDC-6400 COMPUTER. THE B 

APPROXIMATE STORAGE REQUIREMENT FOR THIS PROGRAM IS 14600 (OCTAL). 8 

EXECUTION TIME IS APPROXIMATELY 25 SECONDS PER CASE WITH 180 SURVEY 8 

POINTS, 30 TRAILING SEGMENTS, AND 8 ITERATIONS. NOTE THAT THIS PROGRAM B 

YIELOS A PUNCHEO CARO OECK OUTPUT. B 

9 

INPUT DATA SEQUENCE B 

9 

SPAN, GAMAM, SPEED, ASPECT, NW (4F10. 5,110) B 

SPAN IS WING VORTEX SPAM, FEET B 

GAMAM IS WING CIRCULATION, SQUARE FEET/SE'COND 8 

SPEED IS REMOTE WINO SPEED, FEET/ SECOND B 

ASPECT IS ASPECT RATIO OF THE GEOMETRIC WING BEING REPRESENTED BY B 
THE VORTEX SPAN. VORTEX SPAN IS PI/4 TIMES GEOMETRIC SPAN. 9 

NW IS THE NUMBER OF TRAILING SEGMENTS IN THE WAKE, LESS THAN 50 B 

B 

DELTAX (F10.5) B 

LENGTH OF TRAILING SEGMENTS, FEET. USUALLY TAKEN SPAN/10 9 

B 

XW(i), YW ( 1) (2F10.5) B 

X AND Y COORDINATES OF CENTER OF BOUND VORTEX, USUALLY 0.0, 0.0 B 

X AXIS IS POSITIVE DOWNSTREAM, Y IS POSITIVE UPWARD, 7 TO RIGHT 9 

LOOKING DOWNSTREAM B 

9 

TLMN , TLMX, BELT X (3F10.5) B 

MINIMUM ANO MAXIMUM TAIL LENGTHS, FRACTION OF SPAN, DEFINING 3 

LONGITUDINAL REGION TO 3E SURVEYED, ANO INCREMENT BETWEEN B 

SURVEY POINTS, FRACTION OF SPAN. 9 

B 

THMN , THMX, OELTY (3F10.5) B 

MINIMUM AND MAXIMUM TAIL HEIGHTS, FRACTION OF SPAN, DEFINING B 

VERTICAL REGION TO BE SURVEYED, AND INCREMENT BETWEEN SURVEY B 

POINTS, FRACTION OF SPAM. B 

9 

THSP, OELTZ (2F10.5) B 

SEMISPAM OF TAIL, FRACTION OF SPAN, DEFINING LATERAL REGION TO B 

9E SURYEYEO, AND INCREMENT BETWEEN SURVEY POINTS, FRACTION OF 9 

SPAN. 8 

B 

KK (ID 8 

INTEGER VARIABLE SET EQUAL TO ONE IF SURVEY REGION ABOVE IS 9 

REFERENCED TO WING, ANO TO ANY OTHER VALUE IF REFERENCED TO 8 

SPACE COORDINATES. B 

8 

ADDITIONAL CASES* 9 

REPEAT THE PRECEDING SET OF SEVEN DATA CARDS FOR AS MANY CASES B 

AS DESIREO 9 

!> B 

PUNCHEO OUTPUT RESULTING FROM EACH CASE WILL BE AS FOLLOWS B 


0 

1 

2 

3 

4 

5 

6 
7 
5 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 
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c 



8 

52 

C CARO 1-3. VORTEX SPAN, REMOTE VELOCITY, MING CIRCULATION, ASPECT 

B 

53 

C RATIO, LIFT, DRAG, TOTAL X-VELOCITT AT WING CENTER SPAN, TOTAL 

3 

54 

C Y- 

VELOCITY 

AT WING CENTER, WING GEOMETRIC ANGLE OF ATTACK (4E2U.10) 

B 

55 

C 



3 

56 

C CARD 4 AND 

FOLLOWING CAROS, COORDINATES OF SURVEY POINTS XCI, YCJ, 

3 

57 

C AND ZCJ (SPACE FIXEO) AND TOTAL X, Y, AND Z VELOCITY COMPONENTS 

3 

58 

C AT 

EACH SURVEY POINT. (4E20.:i0J 

3 

59 

C 



3 

60 

C LAST CARO. 

THE NUMBER 10000 IS PUNCHED TO INDICATE THE END OF 

3 

61 

C EACH CASE. 

THIS SPECIAL PUNCHING IS USED BY THE WING-IN-TUNNEL 

3 

62 

C PROGRAM TO 

LOCATE THE ENO OF EACH DATA OBCK. (40X, E20.10> 

B 

63 

C 



B 

64 

C 



8 

65 

1 

FORMAT 

(4F10. 5,11 0) 

B 

66 

2 

FORMAT 

(F10.5) 

B 

67 

3 

FORMAT 

( 18H ITERATION NUMBER ,12 » 

3 

68 

4 

FORMAT 

(10F10.5) 

B 

69 

5 

FORMAT 

(3F10.5) 

B 

70 

6 

FORMAT 

C10F12.6) 

8 

71 

7 

FORMAT 

( 2F10 • 5)i 

B 

72 

8 

FORMAT 

C13H CL/ASPECT = , F8. 5,15 X, 13HCDI/ASPECT = ,F8.5) 

B 

73 

9 

FORMAT 

( 1 3 ,5 F10.5 ) 

B 

74 

10 

FORMAT 

( 2110) 

3 

75 

11 

FORMAT 

(12) 

B 

76 

12 

FORMAT 

(7F15.5) 

B 

77 

4100 

FORMAT 

(74H-N0TE - ALL DISTANCES MEASURED FROW ASSUMEO LIFTING LIN 

B 

78 


IE POSITION AT XHCl) 1 

3 

79 

4150 

FORMAT 

( 18H WAKE COORDINATES ,/, 9X, 2HXW , 13X , 2HYW,13X,2HZW,13X, 

B 

80 


13HDSM1 


B 

81 

4160 

FORMAT 

(4F15.5) 

B 

82 

4170 

FORMAT 

( 1H0, 8HGAM AM = ,F10.4) 

3 

83 

4175 

FORMAT 

(1X,I2,3F10.4,2X,3F10.4,2X ,3F10.4) 

B 

84 

4180 

FORMAT 

( 19H ANGLE OF ATTACK = ,F6:.3,12H RADIANS OR ,F7.3,8H DEGREE 

3 

85 


IS > 


B 

86 

4190 

FORMAT 

(22H ANGLE OF ZERO LIFT = ,F6.3,12H RADIANS OR ,F7.3,8H DEG 

B 

87 


1REES 

) 

3 

88 

42 00 

FORMAT 

C23H TAIL SPAN (ABSOLUTE) =,F9.4,2X,21HTAIL SPAN/WING SPAN 

B 

89 


1=,F9.4> 

B 

90 

4250 

FORMAT 

(1H1»5X,7HSPAN = , F6. 3 , 21X ,18HREMOTE VELOCITY = ,F9.3, 

B 

91 


17 X »14HCIRCULATI0N = ,F9. 3 ,/,6X,15 HASPECT PATIO = ,F6. 3,13X,7HLIFT 

3 

92 


1= ,F9.i 

4,18X,7HDRAG = , F9. 5 ,/ , 6X, 1 3HVX AT WING = ,F10.4,11X, 

B 

93 


113HVY 

AT WING = , F10.4,11X,18HG£OMETRIC ALPHA = ,F6.2,8H DEGREES, 

B 

94 


i/,/,/,1 

BX , 16HWING COORDINATES, 18X, 17HEARTM COOROI NATES ,17X , 

3 

95 


119HVEL0CITY COMPONENTS ) 

B 

96 

42 60 

FORMAT 

(1H , 3F10. 4,5X, 3F10.4,5X,3F1C,4) 

3 

97 

4270 

FORMAT 

(12H TAIL SPAN = , F8.4,4X,'23HTAIL SPAN/WING SPAN * ,F8.4) 

B 

98 

42 80 

FORMAT 

(4E20.10) 

B 

99 

42 81 

FORMAT 

( 4 OX, E20. 1 0) 

3 

100 

42 85 

FORMAT 

(4F10.4) 

B 

101 

42 90 

FORMAT 

( 1H0, 12H2- 0 ALPHA = ,F8.5,12H RADIANS OR ,F7.3,8H DEGREES) 

3 

102 

42 95 

FORMAT 

(1H ,16HIN0UCE0 ALPHA = ,F!8.5,12H RADIANS OR ,F7.3,5H DEG.) 

B 

103 

43 00 

FORMAT 

(1H , 18HGE OMETRIC ALPHA = ,F8.i5,9H RAD. OR ,F7.3,5H DEG.) 

B 

104 

4310 

FORMAT 

( 1 Hi) 

8 

105 

43 2C 

FORMAT 

(1H0,44X,3HX = ,F9.4) 

'9 

196 

43 30 

FORMAT 

(1H0,44X,3HY =, F3. 4) 

B 

107 
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4340 FORMAT (1H , A IX, 1 8HREFERE NCED TOi WING) 9 

4350 FORMAT (1H , 40X, 2 OHREFERE NCED TO TUNNEL) 9 

REAL LIFT 9 

DIMENSION VX(7),VY(7)»VZ(7) 9 

0IMENSI9N VMX(7) , VMY (7) , V MZ( 7) 9 

DIMENSION VCX(7> , VCY(7>, VCZ(7) B 

DIMENSION XW (50) , YW(50>, ZW(50), RW(2,2), DSM(5 0) , VBAR( 2) 9 

DIMENSION ALPHA(7I , BETA (7) B 

RHO = .002378 9 

30 CONTINUE B 

READ (5,1) SPAN,GAMAM,S 3 EEO, ASPECT, NW B 

IF (E0F,;5> 60,31 B 

31 READ (5,2) OELTAX B 

READ (5,7) X W ( 1) , YW( 1) B 

IF (EOF, 5) 60,80 B 

C B 

C COMPUTE INITIAL COORDINATES, WING DIMENSIONS,' TRAILING SEGMENTS B 

80 CONTINUE B 

NW1 = NW 1 B 

ZW(1) = SPAN/2. B 

CHORD = SPAN/ (ASPECT*. 785 398163**2) B 

ALFAA=AS IN (G A MAM* 2./ ( 6. 28 3185 3*CH ORO*SPEE 0) ) B 

XCI = 0.75*CHORO*SQRT(1.- (.785398 16**2)) 9 

XW (2) = XW(1) XCI*COS(ALFAA) B 

YW(2) = YW(1) - XCI*SIN( ALFAA) B 

ZW (2) = ZW(1) B 

XCI = OELTAX ♦ XW(2) 9 

YCJ a YW(2) B 

ZCJ ! ZK(1) B 

DO 90 N= 3 ,NW 9 

ZW (N) a ZCJ 9 

YW(N) a YCJ B 

XW(N) a XCI 9 

XCI = XCI ♦ OELTAX 9 

90 CONTINUE 8 

XW(NWi) a XW(NH) ♦ 1000.0 B 

YW (NW 1) a YCJ B 

ZW(NWl) a ZCJ B 

00 81 1=1, NW 8 

J = I»-1 B 

81 OSM(I) = SQRT( (XW (I)-XW( J) )**2*(YW(I)-YW( J) ) **2+ (ZW (I ) -ZW ( J) ) **2) B 

C B 

C CARRY OUT ITERATIVE SOLUTION 3 

NUMIT = GAMA'M/19. *3. 9 

WRITE (6,4310) B 

DO 100 NUMBER = 1, NUMIT .3 

CALL WKIT ( XW,YW,ZW,OSM, GANAM, SPEED, SPAN, NW,NW1, 9 

1 ALP HAO , ALPH AI , ALFAA , CHORD) 9 

IF ( (NUMIT-NUMBER).GT. 3) GO TO 93 3 

WRITE (6,3) NUMBER 9 

WRITE <5 , 415 0) B 

WRITE (6,4160) (X W(L ) , YW ( L) , ZW (L ) ,OSH (L) , L=l, NW1 ) 8 

CALL LCOMP (XW, YW, Z W,OSM,GA MAM, SPEE 0, SPAN, NW,NW1, LIFT, RHO, B 

1VXWC, VYWC,DRAG) B 

WRITE (6,4170) GAMAM B 

ALPHAO = -ALPHAO 9 
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158 

159 

160 
161 
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ALFAA * -ALFA A 
0£G=ALPHAC*57. 29578 
WRITE (6., 4290) AL PHAQ, DEG 
D£G=ALPHAI*57. 29576 
WRITE (6,4295) ALPHAI,OEG 
OEG = ALFAA* 57. 29 578 
WRITE (5, 4300) ALFAA, OE3 
95 XCI = XW(1) 

00 1000 L = 4 , NW1 
IF (XW(L) .LT.XCI) GO TO 999 
10 00 CONTINUE 
10 0 CONTINUE 
C 

C SET UP COORDINATES FOR VELOCITY SURVEY 
READ (5,5) TLMN,TLMX, CELT X 
READ (5,5) THMN,T HMX ,OELT Y 
REAO (5,7) THSP.OELTZ 
NTL=INT( ( TLMX-TLM N)/OELTX *0. 5 ) 4-1 
NTH=INT( (THNX-THHN)70ELTY+0.5)+l 
NTS=INT( THSP/OELT Z+0.5 ) *1 
COSA=COS (ALFAA) 

SI NA= SIN (-ALFAA) 

WRITE (7,428 0) SP AN, SPEED , GAM AN, A SPECT ,LI FT ,ORAG , VXWC , VYWC, ALFAA 
REAO (5,140) KK 
40 FORMAT (ID 

00 400 1 = 1, NTH 

YC=(THKN*FLOAT(I-D*DELTY)*SPAN 

WRITE (6,4250) SPAN,SPE; D ,GAMAM, ASPECT,LIFT,ORAG , VXWC, VYWC, OEG 
WRITE (6,433 0) YC 
IF (KK.EQ.l) WRITE (6,4340) 

IF (KK.NE.l) WRITE (6,4350) 

00 400 J=i,NTL 

XC= (TLMN +FL0'AT ( J- 1) *OELT X ) *SPAN 
WRITE (6,4320) XC 
IF (KK.EQ.l) WRITE (6,4340) 

IF (KK.NE.l) WRITE (6,4350) 

00 400 K = i ,NTS 
IF (KK.NE.l) GO TO 51 
XCI=XC*COSA*XW(i) -YC*SINA 
YCJ=XC*SINA*YC*COSA*YW(i) 

ZCJ=FLOAT (K-1)*0ELTZ*SPAN 
GO TO 52 

51 CONTINUE 
XCI=XC*XW Cl)t 
YC J=YC*YW (1) 

ZCJ=FLOAT(K-l)*OELTZ*SPAN 

52 CONTINUE 
C 

C COMPUTE VELOCITY COMPONENTS AT SURVEY POINTS 

CALL VCOMP (XCI, Y3J ,ZCJ, DSM,GAMAN, SPAN, SPEED, 

1 VXMOO , VYMOD, VZMOO , VXTOT , V YTOT ,VZFOT ,XW, YW , ZW, NW, . FALSE. ) 

C 

C REFERENCE S«»ACE FIXED COORDINATES TO BOUND VORTEX 
XCI=XCI-XW(1) 

YCJ=YCJ-YW(1) 

WRITE (7,4280) XCI,YCJ,ZCJ, VXTOT, VYTOT.VZTOT 
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216 
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O O 


WRITE (6 i 4260) XC .VC, ZCJ , XCI , YCJ, 2C J, VXTOT , VTTOT , VZTOT B 220 

43 0 CONTINUE B 221 

ZCJaiOOOO. 8 222 

WRITE IT 1 423 1) ZCJ 9 223 

8 224 

READ INPUT DATA FOR NEXT CASE 9 225 

GO TO 30 B 226 

999 CONTINUE B 227 

60 STOP 9 228 

END 8 229 
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SUBROUTINE HKIT ( XH,YH , ZW , OSM, GAN AM , SPEED ,SPAN, N W ,NW1 
i , ALPHAO,ALPHAI,ALFAA ,CHORO) 

C 

C SUBROUTINE TO ITERATE TRAILING WAKE POSITION 
C 

0IMENSI9 N XW (50)»YW(50)»ZW(5Q) , OS M (50 ) »RW (2, 2 » , V BAR (2 ) 

LOGICAL SKP, WTEST 
SKF = .FALSE. 

C 

C MAKE TWO PASSES, FIRST FOR X-Y MOVEMENT, SECOND FOR X-Z MOVEMENT 
DO 20 N = 1,2 

31 NNN = NW 

*0 00 47 M = 1, NNN 

WTEST = .FALSE. 

IF (M.EQ.l) WTEST = .TRUE. 

IF ((M.EQ.l) .AND. SKP) GO TO 47 
C 

C CHOOSE COORDINATES FOR VELOCITY COMPUTATION 

41 XCI = XW(M) 

YCJ = YW(M) 

IF (M.EQ.l) GO TO 42 
ZCJ = ZW (M) 

GO TO 43 

42 ZCJ - 0.0 

43 CONTINUE 
C 

C COMPUTE VELOCITY COMPONENTS AT CHOSEN COORDINATES 

CALL VWKIT (XCI, YCJ, ZCJ, DSM ,GA MAM, SPAN, SPEED, 

1VXMOD, VYMOO, VZMOD , VXTOT, VYTOT ,VZT OT ,XW»YW , ZW,iNW , WTEST ) 

VXM = VXTOT 
VYM = VYTOT 
VZM = VZTOT 

VEL = SQRT (VXTOT**2 ♦ VVTOT**2 ♦ VZTOT**2) 

J = H*1 
C 

C COMPUTE NEW ANGLE OF ATTACK OR SEGMENT ORIENTATION, AND SHIFT 
C TO BE APPLIED TO FOLLOWING SEGMENTS 
IF (M.NE.l) GO TO 45 

ALPHAC=ASIN(-GAMAM*2./(6. 2831 853* CHORD* VE L) ) 

ALPHA I = ATAN ( VYM/ VXM) 

ALFAA = ALPHA3 ♦ ALPHAI 

XSHFT = 0SM(1)*C0S(ALFAA) ♦ XW(1) - XW(2) 

YSHFT = DSM(1)*SIN( ALFAA) ♦ YW(1) - YW(2) 

ZSHFT = 0.0 
GO TO 57 

45 DCWX = VXM/VEL 

XSHFT = OSM( M) *OC MX ♦ XW(M) 

XSHFT = XSHFT - XW(J) 

IF (SKP) GO TO 49 

DCWY = VYM/VEL 

YSHFT = DSM( M)*DCWY ♦ YW ( M) 

YSHFT = YSHFT - YW(J) 

GO TO 57 

49 DCWZ = VZM/VEL 

ZSHFT = OSM(M)*DCWZ ZW(M) 

ZSHFT = ZSHFT - ZW(J) 
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IF (J.EQ.NWD ZSHFT = 0. 

9 

286 

c 


9 

287 

c 

COMPUTE NEW COORDINATES OF TRAILING SEGMENTS DOWNSTREAM OF 

9 

288 

c 

NEWLY ORIENTED SEGMENT 

9 

289 

57 

DO 48 L- J » NW1 

9 

290 


XW <LJ = XW(L) ♦ XSHFT 

9 

291 


IF tSKP) GO TO 59 

9 

292 

59 

YW (L) = YW(L) «• YSHFT 

9 

293 


GO TO 50 

9 

294 

59 

ZW(L> = ZW(L) ♦ ZSHFT 

9 

295 

50 

K = L-l 

9 

296 


OSM(K) = SORT ( (XW (L) -XWC K ) ) **2+< Y W (L) -YW( K) ) **2* C ZW (L ) -ZW {«>) **2> 

9 

297 

48 

CONTINUE 

9 

298 

47 

CONTINUE 

9 

299 

C 


9 

300 

c 

RETURN FOR NEXT PASS 

9 

301 


SKP = .NOT.SKP 

9 

302 

20 

CONTINUE 

9 

303 


RETURN 

9 

304 


END 

9 

305 
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UUU -H CVJ 


SUBROUTINE LCOMP CXW, YW, Z W,DSM,GA MAM,SPEE 0, SPAN, NW, NW1, LIFT ,RHO, 
1 VXTOT,YYTOT,ORAGI 

SUBROUTINE TO COMPUTE LIFT iNO INDUCED ORAG ON WING 

DIMENSION XW<50>, YW( 50 ) , Z W<50 ) ,DSM(50) ,RW (2»2)»VBAR(2),ALPHA(7) , 
1BETA(7> 

real lift 

FORMAT ( 1H0, 6HLIFT =,F13. 4»5X ,6H0RAG =,F1Q.4> 

FORMAT < 1H0, 17HCL/ASPECT RATIO = ,-F10. 4, 5X »17HCD/ ASPECT RATIO = , 
1F10.4) 

FORMAT C 1H0, 23HVX AT WING CENTERLINE = ,Fl 0.5, /,1 HO , 23HVY AT WING 
1ENTERLINE =,F10.5) 

KK = 1 
XCI = XW(KK) 

YCJ = YW ( KK) 

ZCJ = 0. 

CALL VLCOMP < XCI , YCJ, ZCJ, DSM,GAMAM, SPAN, SPEED, 

1YXMOO,VYMOO,VZMOO,VXTOT, VYTOT,VZTOT, XW,YW ,ZW,NW, .FALSE.) 

LIFT = RHO*VXTOT*SPAN*GAMAM 
DRAG = -RHO*VYTOT*SPAN*SAMAM 

CDIAR a ( 3.14159/4 .) **Zt C .5*RHO*< SPEE0**2> MSPAN**2)) 

CLAR = LIFT*COIAR 
CDIAR = ORAG*COIAR 
WRITE (S , 1) LIFT, ORAG 
WRITE (6,2) CLAR, CDIAR 
WRITE C6, 3) VXTOT ,VYTOT 

RETURN 

ENO 
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SUBROUTINE VWKIT <XCI ,YREF, ZREF, 

DSM, GAM AM, SPAN, SPEED 

9 

336 


l,VXMOO,VYMOD ,VZMOO,VXTOT, VYTOT,VZTOT,XW, 

Y W, ZW,NW , WTEST) 

3 

337 

c 




B 

338 

c 

SUBROUTINE TO COMPUTE VELOCITY COMPONENTS 



9 

339 

c 




9 

340 


DIMENSION XM(50), YWC50) , Z W(50 > ,OSM(50) ,RW<2,2), 

V BAR (2) 

B 

341 


LOGICAL WTEST, LTEST 



8 

342 


LTEST = .FALSE. 



B 

343 


GO TO 10 



9 

344 


ENTRY VCOMP 



8 

345 


LTEST = .FALSE. 



B 

346 


GO TO 10 



B 

347 


ENTRY VLCOMP 



9 

348 


LTEST = .TRUE. 



9 

349 

10 

VXM=0.0 



9 

350 


VYM = 0. 0 



9 

351 


VZM s 0.0 



9 

352 


YCJ = YREF 



9 

353 


ZCJ = ZREF 



B 

354 


P = 6.28 31853 



B 

355 

c 




B 

356 

c 

INITIALIZE VARIABLES TO COMMUTE VELOCITY INDUCED BY 

THE SEGMENT 

B 

357 

c 

PAIR UNDER CONSIDERATION 



B 

358 


XHK. s. -XN t ll s 



8 

359 


YHK = YW ( 1) 



9 

360 


ZMK s ZW(1| 



B 

361 


RW12 = CXWK-XCI»**2 ♦ IYW K-YC J>**:2 



B 

362 


RM112 = RW12 ♦ (Z WK-ZC J) * *2 



8 

363 


RH122 = RW12 + (ZWM-ZCJI **2 



9 

364 


RM11 s SQRTf RW112) 



B 

365 


RW12 * S QRT ( RM122 ) 



9 

366 

43 

DO 46 K s 1, NW 



B 

367 


J = KM 



9 

368 


XWJ = XM ( J) 



9 

369 


YW J = YWtJ) 



9 

370 


ZWJ 3 ZW(J) 



B 

371 


RW22 = (XWJ-XCII**2 ♦ (YWJ-YCJI**2 



9 

372 


RW212 3 RW22 *■ (ZWJ-ZCJ»**2 



B 

373 


RW222 3 RW22 ♦ <ZWJ*ZC J) **2 



9 

374 


RW21 = SQRTCRW212) 



9 

375 


RW22 3 S QRT ( RW222 ) 



9 

376 


OSMK 3 OSM(K) 



9 

377 


OSMK2 3 DSMK**2 



9 

378 


H 3 4.*RWH2*oSMK2 - ( RW1 12-RW212 «-OSMK2> 

3* 2 


B 

379 


IF tH.LT.l.E-10) GO TO 44 



8 

380 


VBAR1 3 -GAMAM*<DSMK2-UHU-RH21> **2) *CRW 11+RW21I / CP*RW11*RW21*H) 

B 

381 


GO TO 45 



9 

382 

44 

VBAR1 3 0,0 



8 

383 

45 

H * 4. *RW122*OSMK2“CRW122“RW222*DSMK2) **2 


9 

384 


IF CH.LT . l.E-10) GO TO 4Z 



8 

3 85 


V8AR2 s -GAMAM*(DSM<2- (^W12-RW22) >*2)*CRW12+RW22)/CP*RW12*RW22*H) 

9 

386 


GO TO 48 



B 

387 

4/ 

VBAR2 = 0.0 



9 

388 

40 

CONTINUE 



3 

389 

C 




B 

390 

C 

COMPUTE VELOCITY COMPONENTS INOUCEO BY EACH 

SEGMENT 

PAIR 

9 

391 
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O O <T' CJ O 


VXM S V3AR1* ( (YWK-YCJ>MZWJ-7WK»- (ZWK-7CJ ) MYWJ- YWK) ) 8 

1 -V3AR2*((YWK-YCJ)*(ZWK-ZWJ)-I(-ZWK'-ZCJ)*(YWJ-YWK)) ♦ VXM 9 

VYM = V3AR1M (ZWK-ZCJ) *(XWJ-XWK> - CXWK-XCI ) * (ZWJ- 2WK) ) 8 

1 -VBAR2* ( (-ZW K-ZCJ)M XWJ-XWK) - ( XW K~XC I> ♦ ( ZWK- ZWJ) ) + VYM .9 

IF (LTEST ) GO TO 55 8 

VZM = (V8AR1-VBAR2>*( ( XWK-XCI )♦ (V WJ-YWK) - (YWK-YC’J)* (XWJ-XWK)) «-VZM B 
55 CONTINUE 3 

RW11 = RW 21 8 

RW12 = RW 22 9 

. RWH2 = RW212 B 

RW122 = RW222 9 

XWK = XWJ 9 

YWK = YWJ 8 

ZWK = ZWJ 8 

45 CONTINUE B 

IF (WTEST) GO TO 60 8 

XWK = XWCi) 8 

YWK = YH(1) 8 

ZWK = ZW (1) 8 

HM2 = (XWK-XCI>**2 «■ ( YdK-YCJ)**2 B 

IF (HM2.LT.. 00001) GO TO 60 8 

RM1 = SORT (HM2 > (ZWK-ZCJ ) **2) 8 

RM2 = SORT (HM2 «■ (ZWK*ZC J ) **2) 8 

P = 25.13274 8 

B 

COMPUTE VELOCITY INDUCED BY BOUND VORTEX B 

VXM = GAMAM* (RM14-RM2)* (SPAN**2 - < RH1- RM2) ** 2» *( Y C J-YWK) / ( P* SP AN* 8 

1RM1*RM2*HM2) ♦ VXM B 

VYM = GAMAMMRM1»RM2) * (SPAN**2 -( RM1-RM2) **2) * (XWK-XCI) / (P* SPAN* 8 
2RM1*RM2*HM2I ♦ VYM B 

60 CONTINUE 3 

VXMOD * VXM 9 

VYMOD = VYM 9 

VZMOO = VZM 9 

B 

STORE TOTAL VELOCITIES 8 

VXTOT = VXM ♦ SPEED 8 

VYTOT = VYM B 

VZTOT = VZM 9 

RETURN B 

END 8 
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o n o o o 


APPENOIX C 


PROGRAM TO COMPUTE LINEARIZED MALL INTERFERENCE FACTORS 
FOR TUNNELS OF ARBITRARY CROSS SECTION 


PROGRAM STWKWT (INPUT, OUTPUT, TAPE5=INPUT,TAPE6=OUTPUT) C 0 
C Cl 
C THIS PROGRAM COMPUTES LINEARI ZED MIND TUNNEL WALL INTERFERENCE FACTORS C 2 
C FOR WINO TUNNELS WITH VERTICAL ANO LATERAL PLANES OF SYMMETRY IN THE C 3 
C SPECIAL CASE OF THE MOOEL LOCATED ON THE PLANE OP VERTICAL SYMMETRY. C 4 
C THE MODEL IS A SIMPLE HORSESHOE VORTEX SYSTEM. C 5 
C THE CROSS SECTION OF THE TUNNEL MUST REMAIN CONSTANT OVER THE FULL C 6 
C LENGTH. C 7 
C C 8 
C THIS IS A FORTRAN IV PROGRAM WRITTEN FOR THE CDC 6400 COMPUTER. C 9 
C STORAGE REQUIREMENT FOR THIS PROGRAM IS APPROXIMATELY 46000 (OCTAL) C 10 
C LOCATIONS ON THE CDC 6400. C 11 
C EXECUTION TIME ON THE CDC 6400 IS APPROXIMATELY 95 SECONDS FOR ONE C 12 
C CASE INCLUOING THE MATRIX INVERSION.: C 13 
C C 14 
C INPUT DATA SEQUENCE. C 15 
C C 16 
C TITLE (8 A10) C 17 
C ANY TITLE MAY BE USEO TO ACCOMPANY OUTPUT. C 18 
C C 19 
C MM, NN (212) C 20 
C MM IS THE NUMBER OF COORDINATE PAIRS DEFINING THE COMPLETE CROSS- C 21 
C SECTIONAL SHAPE OF THE TUNNEL. MM CANNOT EXCEED 20. C 22 
C NN IS THE NUMBER OF VORTEX RECTANGLES MAKING UP THE LENGTH OF THE C 23 
C TUNNEL. NN CANNOT EXCEED 25. C 24 
C C 25 
C Y, Z (2F15.5) C 26 
C Y ANO Z ARE THE COORDINATES, IN FEET, OF THE POINTS DEFINING THE C 27 
C SHAPE OF THE TUNNEL, MM CAROS ARE' REQUIRED. C 28 


C THE ORIGIN OF THE COORDINATE SYSTEM IS TAKEN ON THE TUNNEL CENTER C 29 
C LINE WITH X POSITIVE OOWNSTREAM, Y POSITIVE UPWARD, ANO Z POSITIVE C 30 

C TO THE RIGHT LOOKING DOWNSTREAM.. THE FIRST CARD IN THE SEOUENCE IS C 31 

C THE FIRST COORDINATE TO THE RIGHT (POSITIVE Z) OF THE POSITIVE Y C 32 

C AXIS, AND SUBSEQUENT POINTS ARE TAKEN CLOCKWISE AROUND THE TUNNEL. C 33 


C SEGMENT LENGTHS BETWEEN ADJACENT POINTS SHOULD BE EQUAL. C 34 

C C 35 

C OELTAX (F15.5) C 36 

C LENGTH IN FEET OF THE VORTEX RECTANGLES IN THE STREAMWISE C 37 

C DIRECTION. SHOULD BE EQUAL TO THE LENGTH OF SEGMENTS IN THE C 38 

C CROSS -SECTION. C 39 

C C 40 

C SPAN (F15.5) C 41 

C VORTEX SPAN, IN FEET, OF THE WINS • C 42 

C C 43 

C ADDITIONAL CASES C 44 

C REPEAT THE LAST CARD, SPAN (F15.5), FOR AS MANY CASES AS DESIRED. C 45 

C C 46 

1 FORMAT (212) C 47 

2 FORMAT (2F15.5) C 48 

3 FORMAT (F15.5) C 49 

4 FORMAT (4F15.5) C 50 


79 



5 FORMAT (1F10.5) C 

7 FORMAT (3F15. 51 C 

9 FORMAT ( 8 A10 ) C 

210 FORMAT C 1H1, 20X,8A10) C 

211 FORMAT ( 1H0* 53X f 21H1 0 0 E L DATA , / ,/, 25X , 7HSP AN = , C 

1F6.3, 5X, 4HXM = ,F6.,3, 5X, 4HYN = , F6. 3 ,5X, 13HCIRCULATI0N = , C 

2F7.3 ) C 

212 FORMAT ( 1H0, 48X,23HT U M N E L OATA 35X, 9H POINT NO. C 

1,7X,1HY,-9X,1HZ,8X , 14HIENGTH OF SIDE , / » (/»38 X, I 2 »7X, F8. 4, 2X, F6. 4, C 
2 9X » F7 . 4) ) C 

213 FORMAT (1H1,54X,13HR E S U L T S, /,/,5X,llHC00R0lN4TES,5X, C 

ilOHCORR£CTION,6X, 16HT0T4L VELOCIT IES.13X, 25HTUNNEL INDUCED VELOCIT C 
2 IE S » 1 0 X, 2 4HM0 DEL INOUCEO VELOCITI ES ,/ ,4X, 1HX, 5X, 1HY ,5X,1HZ, 7X , C 

33HDEL *6X , 2HVX , 9X, 2HVY,9X, 2HVZ.9X, 3HVXC,8X ,3HVYC, 8X, 3H VZC, 8X,3HVXM, C 
48X,3HVYM,8X,3HVZM ) C 

214 FORMAT ( 1H0, 3F6. 2 ,F3. 3 ,3F 11. 4 ,3F1 1. 4.I3F11 .4) C 

215 FORMAT (/ ,/, 48X,17HSECTI0 N LENGTH = ,F7.4) C 

216 FORMAT (/ ,/,45X,22HCR0S3 SECTIONAL AREA = ,F10.4> C 

INTEGER A,9,C,0,E C 

LOGICAL 0PT1, OPT 2 ' C 

DIMENSION X( 26),Y (20) ,Z<20> ,SINPHI(20) ,C0SPHI (20) ,XCPT(25), G 

lYCPT(ll) ,ZCPT(li) ,SIDE(20) ,CC (100 ,100) , SC 25) , GAN AKt 10 0) , C 

1GAMA ( 25. ID » ZM(2) C 

DIMENSION RC26, 20), HL (25, 20) .H0(2C) ,HYZ(2C) C 

DIMENSION GL ('ll), GO(li) C 

DIMENSION TITLEC8) C 

ID * 26 C 

JD = 25 C 

KD - 26 C 

LO = 11 C 

MD * 100 C 

C C 

C REAO TUNNEL ANO MOOEL DESCRIPTION FROM CARDS C 

34 REAO (5,i9) (TITLE (I), 1=1,8) C 

IF (EOF, 5) 700,35 C 

35 READ (5,1) MM,NN C 

IF ( (MM. GT.20 ) .OR. (NN. GT. 25) ) GO TO 700 C 

N1 = NN ♦ 1 C 

REAO (5, 2) (Y(T) , Z(I) ,1=1 ,MM) C 

READ (5,3) OELTAX C 

C C 

C C 

C COMPUTE THE COORDINATES OF THE TUNNEL. C 

CALL COORD ( X, Y,Z , XCPT , VC PT, ZCPT, S» SINPHI »COSPHl * CELT AX, C 

1SIOE, OPT 1 ,OP'T2 , MM ,NN,LL, K K,N1 ,NK, ID, JD,KO ,LD, AREA) C 

C C 

C GENERATE THE MATRIX OF COEFFICIENTS. C 

CALL MATRIX (X,Y,Z,XCPT, YCPT, ZCPT, SINPHI, COSPHI , SIOE, S,CC , C 

1MM,NN,LL,KK,N1,NK,OPT1,OPT2,R,HL,HD,HYZ,ID,JD,KD,LD,MD) C 

C C 

C C 

C COMPUTE INVERSE OF THE CC MATRIX, STORE RESULT IN CC ARRAY. C 

70 CALL INVR (CC »NK-,MD) C 

C C 

C C 

C REAO MODEL OATA FROM PUNCHED CARDS. C 


51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

35 

36 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 
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75 REAO (5, 31 SPAN 

IF (EOF, 5) 70 Of 80 

80 CONTINUE 

C 

C GENERATE THE RIGHT HAND SIOE OF THE MATRIX EOUATION. 

CALL R!HS CS» AN, XM ,T M , ZM, GAHA M,XCPT,T CPT, ZCPT, SINPHI , 

ICOSPHI.GAMAK, JD,KO,LO,M3,NN,KK> . 

C 

C 

C MULTIPLY RIGHT HANO SIDE BY MATRIX INVERSE, STORE RESULT IN GAMA ARRAY 
M = 0 

DO 150 I = 1,NN 

DO 150 J = 1,KK 

N = M ♦ 1 
XCI = 0.0 
DO 130 K = 1 » NK 

130 XCI = XCI ♦ CC(M, K)*GAMAK (K) 

GAMA(I,J) = XCI 
L a LL ♦ 1 - J 
GAMA(I,L) = -XCI 

IF ((.N0T.0P72).AND. (J.EQ.KK) I GAMA(I,JH> = 0.0 
15 0 CONTINUE 

C 
C 

C WRITE RESULTS OF COMPUTATIONS. 

500 FORMAT (30H1 CALCULATED VORTEX STRENGTHS > 

WRITE (6,500) 

DO 502 J = 1,NN 

WRITE (6,501) (GAMA(J,K), K=i,LL) 

501 FORMAT (/,ilF11.6) 

50 2 CONTINUE 

250 FORMAT (81HOOPT1 = .TRUE. THIS IMPLIES VORTEX SINGULARITY AT TOP A 
1ND BOTTOM CENTER OF TUNNEL ) 

IF (OPT1) WRITE (6,250) 

251 FORMAT (85HQ0PT1 = .FALSE. THIS IMPLIES NO VORTEX SINGULARITY AT T 
10P AND BOTTOM CENTER OF TUNNEL ) 

IF ( • NOT . OPT1) WRITE (6,251) 

252 FORMAT (76HOOPT2 = .TRUE. THIS IMPLIES VORTEX SINGULARITY ON PLANE 
1 OF VERTICAL SYMMETRY ) 

IF (OPT2 ) WRITE (6,252) 

253 FORMAT (80H0OPT2 = .FALSE. THIS IMPLIES NO VORTEX SINGULARITY ON P 
1LANE OF VERTICAL SYMMETRY > 

IF (.NOT. OPT 2) WRITE (6,253) 

AO 00 FORMAT ( 27H1RESUL TANT VORTEX STRENGTHS ) 

AO 02 FORMAT (13H0RING NUMBER ,I2,8X,15HX COORDINATE = , F10. A, 8X ,17HM00 
1EL DISTANCE = , F 1C..4, 8X , 22HMODEL OISTANCE/SP AN = ,Fll.A,(/, 

111F11. 6) ) 

AO OA FORMAT (15HOSECTION NUM3ER ,13,7 ,11F11.6) 

AO 10 WRITE (6,A000) 

AO 15 DO A1A0 L=l,Ni 

AO 20 M=L-1 

AO 25 DO A075 1=1, LL 

AO 30 IF (L-2) A050, A060, AOAO 

AO AO IF (L-Nl) A060, A070, A1AO 

AO 50 GL (I) = GAMA (L,I) 

AO 55 GO TO AO 75 


C 107 
C 108 
C 109 
C 110 
C 111 
C 112 
C 113 
C 114 
C 115 
C 116 
C 117 
C 118 
C 119 
C 120 
C 121 
C 122 
C 123 
C 124 
C 125 
C 126 
C 127 
C 128 
C 129 
C 130 
C 131 
C 132 
C 133 
C 134 
C 135 
C 136 
C 137 
C 138 
C 139 
C 140 
C 1A1 
C 142 
C 1A3 
C 144 
C 145 
C 146 
C 147 
C 148 
C 149 
C 150 
C 151 
C 152 
C 153 
C 154 
C 155 
C 156 
C 157 
C 158 
C 159 
C . 160 
C 161 
C 162 
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40 60 

GL (I) = G AMA'(L ,1) - GAMA ( M, I) 


C 

163 

43 65 

GO TO 4075 


C 

164 

40 70 

GL (I) = -GAMA (M,I ) 


C 

165 

43 75 

CONTINUE 


C 

166 

40 77 

XOR = X(L)-XM 


C 

167 

40 78 

XCI e XO R/SP’AN 


c 

168 

43 80 

WRITE (5,4002) L, X (L > , XDR , XC I , (GL (I ) , 1=1 , LL) 


c 

169 

4100 

IF (L-Nl) 4110, 4140, 4140 


c 

170 

4110 

DO 4125 1=2, LL 


c 

171 

4115 

J=I-i 


c 

172 

4120 

GL(J) = G AMA (L , J) - GAMA (L,I) 


c 

173 

4125 

CONTINUE 


c 

174 

4130 

MMM = LL - 1 


c 

175 

4135 

WRITE (5 , 4004) L, (GL ( J) , J = l, MMM) 


c 

176 

4140 

CONTINUE 


c 

177 


WRITE (6,210) TITLE 


c 

178 


WRITE (5,212) (I,Y(I),Z(I) , SIOE( I ) , 1=1, MM) 


c 

179 


WRITE (6,215) OELTAX 


c 

180 


WRITE (6,216) AREA 


c 

181 


WRITE (6,211) SPAN,XM,YM, GAMAM 


c 

182 

C 



c 

183 

c 



c 

184 

c 



c 

185 

C NOW 

BEGIN SURVEY OF TUNNEL FLOW FIELD. 


c 

186 

C PERFORM SURVEY IN THE PUNE OF THE MODEL.* SURVEY FROM APPROXIMATE 

c 

187 

C GEOMETRIC WINGTIP TO CENTERLINE OF TUNNEL WITH FIXED X COORDINATE, 

c 

188 

C THEN 

1 SURVEY ALONG CENTERLINE OF TUNNEL DOWNSTREAM FROM BOUND 

VORTEX. 

c 

189 

C SURVEY INCRE MENT IN 30TH DIRECTIONS IS (VORTEX SPAN)/20 


c 

190 

C SURVEY 9EGINS AT SOUND VORTEX AND CONTINUES FOR THREE VORTEX 

SPANS 

c 

191 

C OOWNSTREAM OF THE SOUND VORTEX. 


c 

192 

C 



c 

193 


WRITE (6,213) 


c 

194 


OTP = SPAN/20.0 


c 

195 

C SET 

XTP, YTP, ZTP TO INITIAL SURVEY COORDINATES. 


c 

196 


XTP = XM 


c 

197 


YTP = YN 


c 

198 


ZTP = S»ANM3./20. 


c 

199 

60 0 

CONTINUE 


c 

200 


CALL SURVEY ( XTP, YTP, ZTP, X , Y , Z,XM , YM, ZM, SINPHI,COSPHl ,S, 


c 

201 

IGA MA, SIDE, 0»T1, SPAN, GAMAM ,VXC,VYC ,VZC,VXT , VYT ,VZT ,VXM , VYM 

»VZM, 

c 

202 

1LL,MM,NN,N1,R,HL,H0,HYZ,I0,J0,KD,LD > 


c 

203 


XOR = XTP-XM 


c 

204 


DEL = VYC*AREA/SPAN/GAMAM/2. 


c 

205 

C 



c 

206 

C V*T 

ARE TOTAL VELOCITY COMPONENTS (SUM OF V*C AND V*M>. 


c 

207 

C V*C 

ARE VELOCITY COMPONENTS INOUCED 3Y TUNNEL WALLS. 


c 

208 

C V*M 

ARE VELOCITY COMPONENTS INDUCED BY MODEL. 


c 

209 

C XOR 

IS X COORDINATE OF SURVEY POINT RELATIVE TO BOUNO VORTEX. 


c 

210 


WRITE (6,214) XOR,YTP,ZTP,OEL,VXT,VYT,VZT,VXC,VYC,VZC,VXM 

, VYM , VZM 

c 

211 


IF (ZTP. GT.0.0) GO TO 601 


c 

212 


XTP s XTP t OTP 


c 

213 


ZTP = 0. 0 


c 

214 


IF (XTP.LE.XM«-3.0*SPAN) GO TO 600 


c 

215 

C 



c 

216 

C READ DATA FOR NEXT MODEL FROM PUNCHE3 CAROS. 


c 

217 


GO TO 75 


c 

218 
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c c 219 

601 ZTP = 7TP - OTP C 220 

GO TO 60 0 C 221 

63 3 CONTINUE C 222 

700 STOP C 223 

ENO C 224 
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SUBROUTINE COORD (X,Y,Z,XCPT, YCPT ,ZCPT,S, SINPHI, COSPHI,DELTAX, 

1 SIDE i OPT 1 , OP"T2*MM , NN,LL,KK,N1 ,NK, ID,JD,KD ,LO, AREA) 

C 

C THIS IS A SUBROUTINE TO COMPUTE THE TUNNEL COORDINATES. 

C 

LOGICAL 0PT1.0PT2 

OIMENSION X(IO),Y (KO> ,Z(KO) ,XCPT( JD), YCPT (LD> , ZC PT( LO) ,S( JD), 
lSINPHI(KO) ,COSPHI <K0> ,SIOE(KD) 

C 

C COMPUTE VORTEX RING X-COOROIN ATES. 

XCI = 0. 0 

00 20 I=i»NN 
XCI) = XCI 

20 XCI = XCI ♦ OELTAX 

X CNi) = 1 GOJ • 0 ♦ X(NN) 

C 

C TEST TUNNEL SHARE COORDINATES ANO OET ERMINE TOTAL NUMBER OF 
C UNKNOWNS CNK). 

OPT1 a Z (MM) .EQ.Q.O 

1 = MM/A 

J s (MM/A) ♦ 1 

0PT2 = ( Y (I) .EQ.O.C) .OR.,(Y (J) .EQ. 0.07 
IF (.NOT.OPT1) GO TO 10 
LL = MM/2 
KK a MM/A 
GO TO 1A 

10 IF (. NOT .OPT 2) GO TO 12 
KK * MM/ A ♦ 1 
GO TO 13 

12 KK a MM/ A 

13 LL = MM/ 2 ♦ 1 

1A CONTINUE' 

NL = NN * LL 

NM a NN* MM 
NK — NN*KK 

IF (NK.LE.100) GO TO 17 
C 

C IF NK IS GREATER! THAN 100, TERMINATE EXECUTION. 

WRITE (5,15) NK 

15 FORMAT ( 1H0, 25M0I MENSIONS EXCEEDED, NK =,I3,16H REDUCE MM OR NN 
STOP 
C 

C GENERATE VORTEX RECTANGLE PARAMETERS. 

17 00 21 I - 1, NN 

21 S(I) = X(I»1) - X(I) 

DO 23 1=2, MM 

22 SIOE(I) = SORT ( (Y (I) - Y(I-1»)**2 ♦ (Z(I) - Z(I-i))**2) 
SINPHI(I) a ( (Y(I)-Y(I-l) )/(SI0E(I))) 

23 COSPHI(I) = ( (Z(I)-Z(I-l) >/(SIDE(I))) 

SIOE(i) = SORT ( (Y (1) - X ( MM) ) **2 ♦ (Z(l* - Z(MM) )»*2) 

SINPHI(l) a ( (Y(l)-Y(MM) ) / (SIDEd )) ) 

COSPHI(l) a ((Z(l)-Z(MM) )/(SIDE(l))) 

C 

C GENERATE CONTROL POINT LOCATIONS. 

DO 2A I = 2, LL 

YCPT(I) = CY(I>*V(I-l)>/(2.) 


C 225 
C 226 
C 227 
C 228 
C 229 
C 230 
C 231 
C 232 
C 233 
C 23A 
C 235 
C 236 
C 237 
C 238 
C 239 
C 2 AO 
C 2A1 
C 2A2 
C 2 A3 
C 2AA 
C 2A5 
C 2A6 
C 2A7 
C 2A8 
C 2A9 
C 250 
C 251 
C 252 
C 253 
C 25A 
C 255 
C 256 
C 257 
C 258 
C 259 
C 260 
C 261 
C 262 
C 263 
) C 26A 
C 265 
C 266 
C 267 
C 268 
C 269 
C 270 
C 271 
C 272 
C 273 
C 27A 
C 275 
C 276 
C . 277 
C 278 - 
C 279 
C 280 
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OOO U o o 


24 ZCPT(I) = <ZCI>+Z(I-l>>/<2.> C 281 

ZCPT(l) = CZ tl>*Z (Hi) )/(2.) C 282 

YCPTUI = (Y(i»+Y<MMI)/C2.» C 283 

MMM = NV - 1 C 284 

DO 25 I * 1, HMM C 285 

25 XCPT til = (X(WI ♦ XUn/<2.) C 286 

XCPT(NN) = X <NN) ♦ DELTAX/2.0 C 287 

C 288 

GENERATE TUNNEL CROSS SECTIONAL AREA. C 289 

AREA =0.0 C 290 

J = MM C 291 

00 30 I = 1, MM C 292 

AREA = AREA * A3S IYCI)-Y< J) ) ♦ ABS ( Z (II +Z C J ) » C 293 

J = I C 294 

AREA = AREA/2. C 295 

C 296 

RETURN TO CALLING PROGRAM. C 297 

C 298 

RETURN C 299 

END C 300 
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SUBROUTINE MATRIX (X , Y ,Z , XCPT .YC* T, ZCPT ,S INPHI, C OSPHI , SIDE, 

S, CC , 

c 

301 


1MM,NN,LL,KK, N1,N< ,0?T1 ,3P T2, R, HL, HO, HYZ,I D, JD ,KD ,LO ,MD> 


c 

302 

c 


C 

303 

c 

THIS IS A SUBROUTINE TO GENERATE THE MATRIX OF COEFFICIENTS FOR 

-1 

X 

in 

c 

304 

c 

SPECIAL CASE OF VERTICAL SYMMETRY. 


c 

305 

c 



c 

306 


LOGICAL OPT1 , OPT2 


c 

307 


INTEGER A,B,C,D,E 


c 

308 


DIMENSION X( ID) ,f (KO> , Z( K 0) , SINPH I ( KD) ,COSPHI (KD) ,XCPT(JO) , 


c 

309 


lYCPT(LO) ,ZCPT(LO> ,R(I0 ,<0 1 ,SIOE(K D) ,CC CMD ,M0) , HL (ID»KQ) , HO(KD) , 

c 

310 


1S( J0> ,HYZ (KO) 


c 

311 


P = 25.13274 


c 

312 

c 

CYCLE THROUGH CONTROL POINTS. 


c 

313 


M = 0 


c 

314 


00 50 1=1, NN 


c 

315 


DO 49 J = 1, KK 


c 

316 


M = M ♦ 1 


c 

317 

c 



c 

316 

c 

SELECT VARIABLES FOR THIS CONTROL POINT. 


c 

319 


SIN J = S INPHI ( J) 


c 

320 


COSJ = C OSPHI (J) 


c 

321 


xci = xcptci) 


c 

322 


YCJ = YCPT(J) 


c 

323 


ZCJ * ZCPT Ml 


c 

324 

c 



c 

325 

c 



c 

326 

c 

GENERATE COORDINATES OF VORTEX RECTANGLES RELATIVE TO PRESENT CONTROL 

c 

327 

c 

POINT. 1 


c 

328 


DO 26 JJ=1,MM 


c 

329 


HO(JJ) = SORT ( (YCJ-Y (JJ)1 **2 ♦ (ZCJ-7 (JJ> >**2) 


c 

330 


HYZ(JJ) = SQRT (((ZC J-Z( JJ) > *SINPHI( JJ> - (YCJ-Y <JJ))*COSPHI <JJ> ) **2) 

c 

331 


00 26 11*1, Ni 


c 

332 


R(II, JJ) =SQR'T ( (XCI-X (II) ) ** 2 + (YCJ-Y ( J J) ) *♦2*1 ZCJ- Z ( JJ) > 2) 


c 

333 

25 

HL(II,JJ)=SQRT((X (IIJ-XCI)**2 ♦ H YZ (Ji J) ** 2) 


c 

334 

C 



c 

335 

c 

CYCLE THROUGH VO'RTEX UNKNOWNS. 


c 

336 


N = 0 


c 

337 


DO 46 K= 1 , NN 


c 

338 


DO 47 L=i,KK 


c 

339 


N = N ♦ 1 


c 

340 

c 



c 

341 

c 



c 

342 

c 

SELECT VARIABLES' FOR THIS PARTICULAR RECTANGLE OR RECTANGLES. 


c 

343 


0 = L 


c 

344 


E = K+l 


c 

345 


MNIMIZ = 0 


c 

346 

101 IF (OPT1 ) GO TO 15 


c 

347 


A = Q-l 


c 

348 


C = 2*LL-B 


c 

349 


0 = C-l 


c 

350 


IF (B-l) 50,29,27 


c 

351 

27 

IF (LL-B) 50,29,28 


c 

352 

15 

IF (B-l) 50,16,17 


c 

353 

IT 

IF (LL-B) 50,19,11 


c 

354 

11 

A = 9-1 


c 

355 


C = MM-A 


c 

356 
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r 



D = MM-B 



C 

357 


GO TO 28 



C 

358 

18 

A = MM 



C 

359 


C = MM 



C 

360 


0 = MM-t 



C 

361 


GO TO 28 



C 

362 

19 

A = LL-1 



C 

363 


C = LL*1 



C 

364 


0 = LL 



C 

365 


GO TO 28 



C 

366 

28 

RKA = R( K , A) 



C 

367 


RKC = Rl K,C) 



C 

368 


REA = R(E,A) 



C 

369 


REC = RfE ,C) 



C 

370 


HLKC * HL (K,C> 



C 

371 


HLEC = HL <E,C> 



C 

372 


HO A = HD ( A) 



c 

373 


HOC = HO ( C) 



c 

371, 


VA = Y (A ) 



c 

375 


ZA = Z (A ) 



c 

376 


ZC = 2 €C > 



c 

377 


HYZA = HYZ(A) 



c 

378 


HY ZC * HYZ(C) 



c 

379 

29 

SI NL = SINPHICB) 



c 

380 


COSL = COSPHI(B) 



c 

381 


RKB = R( K , 8) ! 



c 

382 


RKD * R<K,D) 



c 

383 


REB = R( E » B) 



c 

384 


RED = RfE,0» 



c 

385 


HL KB = HL <K,B> 



c 

386 


HLEB = HL(E,B) 



c 

387 


HDB = HD ( B) 



c 

388 


HOO = HD ( D) 



c 

389 


SIDEB = SIDE (B) 



c 

390 


OK=S ( K) 



c 

391 


Y9 = Y«) 



c 

392 


ZB = Z(3 ) 



c 

393 


ZO = ZCO) 



c 

394 


XK = X(K) 



c 

395 


XE = X<£ ) 



c 

396 


HYZB = HYZ(B) 



c 

397 


HYZO = HYZ(O) 



c 

398 

C 




c 

399 

C 




c 

400 

c 

COMPUTE VELOCITY COMPONENTS 

INOUCEO BY RECTANGLE OR RECTANGLES, 

c 

401 

c 

TAKE ANY SPECIAL CASES INTO 

ACCOUNT. 

c 

402 


IF (COSJ.EO. 0. 00C0C) GO 

TO 

35 

c 

403 


IF C9-1) 50,16,31 



c 

404 

31 

IF (LL-B) 50,16,32 



c 

405 

IS 

IF (. NOT. OPT1) GO TO 33 



c 

406 

32 

IF CCOSL.EQ. 0. 0090C) GO 

TO 

62 

c 

407 


VY= (COSL/ (P*SI0E3>*<-( (RKA«-RKB)MSIDEB**2-(RKA-RKB)**2)/( ( 

c 

408 


2HLKB**2> *RKA*RKB) ♦ (RKD ♦RKC) MSI DEB* *2 - (RKC-RK D> **2)/ < CHLKC»*2> 

c 

409 


2*RKC*RKD) >MXK-XCI> ♦ ( < REA+RE3) * (SIDES** 2 -( RE A -REB) **2> / C l 

c 

410 


2HLEB**2> *REA*REB> ♦ (REC* REO) * (SI DE B**2 - (REC-RE 0) **2>/ ( (HLEC**2> 

c 

411 


2*REC*RE0) >*(XE-XCII) ♦ 

1. 

/ ( P*QK )* ( t (RKB *REB ) * { 0K**2 -< 

c 

412 
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2RKB-REB) **2)/ ( (HD8**2> *RK B*REB>> * (ZB- ZCJ) - ( (RKD*REO>* ( DK**2-( RKD- 
2RE0>**2) / ( (HOD** 2 ) *RKO*RE 0») * (ZD- ZCJ) ♦ ( (RKC+REC) *(DK**2-(RKC-REC) 
2**2)/<(H0C**2)*RKC*REC>) * C ZC- ZCJ* -C (RKA*REA) • (DK**2- (RK A- RE A) **2) 
2/ ( (HD A** 2) *RK A*R£ A) ) * ( ZA- ZCJ) )) 

GO TO 36 

62 VY = (l./(P*DK>*(( (RKB+RE8)*(OK**2 -( 

2RKB-REB) **2) / ( (HO 8**2) *RK B*RE 9) I * (ZB-2CJ) - ( (RKD* RED) * (DK**2-( RKO- 
2 RED)** 2)/ ( (H 00**2 ) *RKO*R£ 0) ) * (ZD- ZC J) ♦ ( (R KC+REC) * (0K**2- ( RKC- REC) 
2**2)/ ( (HOC** 2) *RKC*REC) ) * (ZC- ZCJ) -( CRKA+REA)* (OK* *2- (RKA- RE A) **2) 
2/( (H0A**2) *R<A*RE A))*(ZA-ZCJ) )) 

GO TO 36s 

33 IF (COSL.EQ. 0.000 00) GO TO 63 

VY = (COSL/('P*SIOE8) *(-( ( RK0*RK9) *(SI0EB**2-( RKQ-RKB) **2)/( ( 

2HL KB**2) *RKD*RK3) ) *( XK-XC I) * ({R£D*R£B> * (SIDEB**2 -(REO-REB) ** 

2 2) / ( ( HLE B**2 ) *REO *REB) )* (XE-XCI)) ♦ 1./ (P*DK> *( ( (RK3*REB) *(0K** 
22- (RKfi-REB)**2)/( (H0B**2) *RKB*REB > ) * ( 7B-Z C J)- ( (RKO*REO> * (OK** 2- 

2 (RKO- RED ) **2 ) / ( (H 00**2) * RKO* RED) ) * ( ZD-2CJ ) ) ) 

GO TO 36 

63 VY = (l./(P*DK>*(( (RKB*REB)*(DK** 

22- (RK8-REB)**2)/( (H08**2) *RKB*RE9 ) ) * C ZB-Z C J) - { (RKD*RED) * (DK**2- 
2 (RKO-REO ) **2 ) / ( (H 00**2) * RKO*REO) ) *(ZO-ZCJ))) 

GO TO 36 

35 VY - 0.30000 

36 IF (SINJ.EQ. 0.00000) GO TO 42 

IF (8-1) 50,55,38 

38 IF (LL-B) 53 ,55,3 9 

55 IF ( • NOT • OPT 1) GO TO 40 

39 IF (SINL.EQ. 0.003UG) GO TO 64 

VZ = (SI Nt_/ ( P*SIOEB) * ( (( RKA* RKB) * (SIDES** 2 -( RKA -RKB) **2) / ( ( 
3HLKB**2) *RKA*RK9) - ( RKC* RKO) * (SI OEB* *2 - (RKC-RK 0) **2 ) / ( ( HL KC** 2) 
3*RKC*RK0) )*( XK-XC I) * ((REC*REO>* (SI0E9**2 -(REC-REO) **2>/( ( 
3HLEC**2) *R£C*RED) - (REA* REB) * (SI DEB**2 - (REA-RE B )**2)/ ( (HLEB**2) 
3*REA*RE3) )*(XE-XCI)) * l. / (P* OK) * ( ( (RKA+REA) * (DK **2 -(RKA 
3-REA) **2) / ( ( H9A** 2)*RKA*REA) - (RKC*REC) * (DK**2 - (RKC-REC)**2) / ( ( 
3HDC** 2) * RKC* REC) ) * (YA-YO J ) * ( ( RKO* RED) * ( 0K**2 - ( RKD-RED) **2) / 

3 ( (H00**2 ) *RKO*RED) - (RKB *REB ) *(0 K**2 - (R KB-REB) **2)/ ( (HDB**2) * 
3RKB*REB) ) MY9-VCJ))) 

GO TO 43 

64 VZ = (t . / (P* DK) * ( ( (RKA+REA) * (OK**2 - (RKA 

3-REA) **2) / (f HOA** 2)* RKA* RE A) - (RKC*REC)* (0K**2 - (RKC-REC)**2)/ ( ( 
3H0C** 2) *RKC*REC) ) MYA-Y3J) ♦ ( ( RKO*REO) * ( DK**2 - (RKD-RED) ** 2) / 

3 ( (HDQ**2) *RKO*R£D ) - ( RKB *REB) *(D K**2 -IRKB-RE9) **2)/ { (HDB**2) * 
3RK8*REB) ) *(Y9-YCJ))> 

GO TO 43 

40 VZ a (1,/(P*0K)*( (( f RK9*RE0)*(0K**2 - ( RKO- RED) **2 ) / ( (HOO** 2) *RKD* 

3 RED) -(RK9 * REB) *(OK**2- (RKB-REB )**2)/ {(HD9**2) *RKB*REB) )* 

3 (Y B-YC J) ) ) 



GO 

TO 43 

42 

VZ 

= 0.00000 

43 

IF 

(MNIMIZ) 50,105,106 

10 5 

B = 

LL*1- 8 


MNIMIZ = 1 

c 

C STORE NORMAL VELOCITY IN CC ARRAY, ACCOUNT FOR VERTICAL SYMMETRY. 
CC1 = VY *COSJ - V Z*SIN J 
GO TO 131 


C 413 
C 414 
C 415 
C 416 
C 417 
C 416 
C 419 
C 420 
C 421 
C 422 
C 423 
C 424 
C 425 
C 426 
C 427 
C 428 
C 429 
C 430 
C 431 
C 432 
C 433 
C 434 
C 435 
C 436 
C 437 
C 438 
C 439 
C 440 
C 441 
C 442 
C 443, 
C 444 
C 445 
C 446 
C 447 
C 448 
C 449 
C 450 
C 451 
C 452 
C 453 
C 454 
C 455 
C 456 
C 457 
C 458 
C 459 
C 460 
C 461 
C 462 
C 463 
C 464 
C 465 
C 466 
C 467 
C 468 
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10 6 CCCM,N> = CC1 - V Y*C0S J *■ V Z*SINJ C 469 

C C 470 

47 CONTINUE C 471 

48 CONTINUE C 472 

49 CONTINUE C 473 

50 CONTINUE C 474 

C 475 

THE MATRIX IS COMPLETE, RETURN TO CALLING PROGRAM, C 476 

C 477 

RETURN C 476 

ENO C 479 
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SUBROUTINE INVR(A ,N, ISIZE > 

THIS IS A SUBROUTINE TO INVERT THE MATRIX A. 

THE INPUT MATRIX A IS OESTROYEO AND REPLACED BY TTS INVERSE. 
A IS ASSUHED TO CONTAIN N ROWS AND COLUMNS OF OAT A. 1 ' 

A IS ASSUMED TO BE DIMENSIONED ISIZE BY ISIZE. 


DIMENSION IPIVOT( 103) , AC I SIZE ,ISI ZE > ,'INOE X (10:0, 2 > .PIVOT <1Q0> 
EQUIVALENCE (IROW.JROW) , ( ICOLUM, J COLUH) , ( AMAX ,T, SWAP) 


15 OO 20 J=1,N 
20 IPI VOT (J ) =0 
30 OO 550 1 = 1, N! 

SEARCH FOR PIVOT ELEMENT 

40 AMAX= 0.0 
45 DO 105 J = l.N 

50 IF (IPIVOTU)-l) 60, 105, 60 
60 DO Iti 0 K = i»N 

70 IF (IPIVOT(K)-l) 80, 100, 740 
80 IF (A3S(AMAX) - ASS ( A ( J, K) ) ) 85,100,100 
85 IROW=J 
90 ICOLUM=K 
95 AMAX=A(J,K) 

100 CONTINUE 
135 CONTINUE 

110 IPI VOT (ICOLUM) =IP I VOT C ICOLUM) +1 

INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

130 IF (IROW-ICOLUM) 143, 250, 140 

140 CONTINUE 

150 DO 200 L = 1,N 

160 SWAP=A(IROW,D 

170 A(IROW,L)=A( ICOLUM.L) 

200 A(ICOLUM,L)=SWAP 

260 INDEX (I, l)=IROW 

270 INDEX (I, 2)=IC0LUN 

310 PIVOT (I) =A(ICOLUM, ICOLUM) 

DIVIDE PIVOT ROW BY PIVOT ELEMENT 

330 A (ICOLUM, ICOLUM) = 1.0 
340 OO 350 L=1,N 

350 A (ICOLUM, L > = A < ICOLUM, L)X. PIVOT (I> 

REDUCE NON-PIVOT ROWS 

380 OO 550 Ll=i,N 

390 IF (Li-ICOLUH) 400, 550, 400 

400 T=A(L1, ICOLUM) 

420 A (|_1, ICOLUM) =0.0 
430 DO 450 L=1,N 


C 480 
C 481 
C 432 
C 483 
C 484 
C 485 
C 486 
C 487 
C 488 
C 489 
C 490 
C 491 
C 492 
C 493 
C 494 
C 495 
C 496 
C 497 
C 498 
C 499 
C 500 
C 5 01 
C 502 i 
C 503 ll 
C 534 ? 

C 505 
C 506 
C 507 
C 503 
C 509 
C 510 
C 511 
C 512 
C 513 
C 514 
C 515 
C 516 
C 517 
C 518 
C 519 
C 520 
C 521 
C 522 
C 523 
C 524 
C 525 
C 526 
C 527 
C 528 
C 529 
C 530 
C 531 
C 532 
C 533 
C 534 
C 535 
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450 AtLl»L>=AlLl»L>-A CICOLUi,U*T 
550 CONTINUE 


INTERCHANGE COLUMNS 

600 00 710 1 = 1, N 
610 l=N+l-I 

620 IF (IN0EX(L, 1) _ IN0EX(L,2) ) 630, 710, 630 
630 JR0W=IN9£X(L ,1) 

640 JC0LUM=IN0EX (L,2) 

650 00 705 K=1,N 
660 SWAP= A CK , JRO'W) 

670 A(K, JROW»=ACK, JCOLUM) 

700 A ( K, JCOL UM) = SWAP 
705 CONTINUE 
710 CONTINUE 
740 RETURN 
ENO 


C 536 
C 537 
C 538 
C 539 
C 540 
C 541 
C 542 
C 543 
C 544 
C 545 
C 546 
C 547 
C 548 
C 549 
C 550 
C 551 
C 552 
C 553 
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SUBROUTINE RHS fSP AN, XM, 1 M , ZM, GAHA M, XCPT, YCPT, ZCP T ,SINPHI , C 

ICOSPHI.GAMAK, J0,KQ,L9, MD,NN,KK> C 

C C 

C THIS IS A SUBROUTINE TO COMPUTE THE RIGHT HAND SIDE OF THE C 

C MATRIX EQUATION FOR THE STRAIGHT WAKE IN WIND TUNNEL PROGRAM. C 

C C 

C C 

DIMENSION XCPT{JO),VCPT<LO),ZCPT(LO>,SINPHI(KO),COSPHICKD>, C 

1ZM(2> ,GAMAK(MO)l C 

c c 

C GENERATE MODEL COORDINATES FOR USE IN GENERATING THE GAMAK MATRIX AND C 
C FOR LATER USE IN THE SURVEY SUBROUTINE. C 

GAMAM = 1.0 C 

I = NN/2 ♦ 1 C 

XM = XCPT (I) C 

YM s 0.0 C 

ZM(t) = SPAN/2.: C 

ZM(2> = -ZM<i> C 

ZM1 = ZN(1) C 

ZM2 = ZM<2> C 

C C 

C GENERATE THE RIGHT HANO SIDE OF THE MATRIX EQUATION. C 

P = 25.13274 C 

C C 

C CYCLE THROUGH CONTROL POINTS. C 

M = 0 C 

DO SO I = 1, NN C 

DO 59 J = 1,KK C 

M = M ♦ 1 C 

c c 

c SELECT VARIABLES FOR THIS CONTROL POINT. C 

SINJ = SINPHI(J) C 

COSJ = COSPHI(J) C 

XCI = XCPT (I) C 

YCJ = YCPT(J) C 

ZCJ = ZCPT(J>I C 

C C 

C COMPUTE VELOCITY INDUCED AT CONTROL POINT BY MODEL. C 

RM1 = SQRT((XM-XCI>**2 (YM - Y5J>**2 ♦ (ZMCl) - ZCJ>**2> C 

RM2 = SQRTC(XM-XCI>**2 ► (YM - YCJ>**2 * (7M(2)-ZCJ)**2> C 

HM1 = SORT (CYCJ-YMI**2 * (ZCJ - ZM(1>>**2) C 

HM2 = SORT <( ! YCJ - YM>**2 ♦ (XCI-XM)**2> C 

HM3 = SQRT(( YCJ-YM)**2 » ( ZCJ-ZMC 2>>**2> C 

IF (COSJ. EQ.O. 00000) GO TO 51 C 

VYM = GAM AM* ( (RM1 ♦PM 2) *(SPAN**2 - (RM1-RM2) **2 > *( XM-XCI) / ( P*S» AN* C 

2RM1*RM2* (HM2**2>) *2./P*( ( 1.* (XCI- XM)/ (RM1 ) )* (ZCJ-ZMi) / (HM1**2) ♦ C 
2 (1. * ( XCI - XM) / (RM2))*(ZM2-ZCJ»/(HM3**2 )>) C 

GO TO 52 C 

51 VYM=G. 00000 C 

52 IF (SINJ. EQ. 0. 00000) GO TO 53 C 

VZM = GAMAM* ( (YCJ-YM) *2.I/P)*( (1, ; ♦( XOI-XM )/RM2) / (HM3**2) - (1. M C 

3XCI-XM)/ CRM1) >/(HMl**2>> C 

GO TO 54 C 

53 VZM = 0. 00000 C 

C C 

C STORE NORMAL VELOCITY COMPONENT IN GAMAK ARRAY. C 


554 

555 

556 

557 

558 

559 

560 

561 

562 

563 

564 

565 

566 

567 

568 

569 

570 

571 

572 

573 

574 

575 

576 

577 
576 

579 

580 

531 

532 

583 

584 

585 

586 
537 
5 88 

589 

590 

591 

592 

593 

594 

595 

596 

597 

598 

599 

600 
601 
602 

603 

604 

605 

606 

607 

608 
609 
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S'* GAMAK <H) = VZM*SINJ - VTH*COSJ 

59 CONTINUE 

60 CONTINUE 

RIGHT HAND SIOE IS COMPLETE, RETURN TO CALLING PROGRAM. 
RETURN 
END 


C 610 
C 611 
C 612 
C 613 
C 614 
C 615 
C 616 
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SUBROUTINE SURVEY (XTP ,r TP, ZTP,X,:Y, 2, XM, Y M, ZM ,SI NPHI, COSPHI ,S, 
1 GAMA , SID E » OPT i,S° AN, GA MAM , VXC , VYC , VZC , VXT , VYT , VI T ,VXM , VYM , V ZM, 
1LL,MM,NN,N1,R,HL, HO, HY Z, I 0 , JO ,KO, LO ) 

C 

C THIS IS A SUBROUTINE TO COMMUTE VELOCITY COMPONENTS AT COORDINATES 
C XTP, YTP , ZTP • 

C 

LOGICAL OPT1 
INTEGER A , 3,C , 0,E 

DIMENSION XtID) ,Y (KD) ,Z(KO> ,SINPHI(KO) ,COSPHI(KD) , 

1 R(ID,<0) ,SIDE(KO) ,HL(IO ,KO> ,HD(KD> ,S(JD) , 

iGAMAt JO,LO), ZMC2) ,HM< 3) , H YZC K D) , RM (2) 

C 

C OEFINE POSITION OF MODEL ANO VORTEX RECTANGLES RELATIVE TO SURVEY 
C POINT, 

ZM1 = ZH <11 
ZM2 = ZM C 21 

601 RM(1» = SORT ( (XM-XTP) **2 + (YM - YTP) **2 ♦ CZMC1) - ZTP>**2) 

RM (2) = SORT ( (XM-XTP>**2 ♦ { YM - YTP)**2 ♦ (ZM(2 ) -ZTP > **2> 

HM1 s SORT ( ( YTP-Y M) **2 t- (ZTP - ZM(1))**2) 

HM2 = SORT ((YTP - YM) **2 ♦ (XTP-XM)**2) 

HM3 = SORT (('YTP-Y M)**2 ♦ (ZTP-ZM( 2) )**2) 

DO 127 J = 1 , MM 

HD ( J) = SQRT((YTP-Y(J))**2 + (ZTP - Z(J))**2) 

HYZ(J) = SQRT(((ZTP-Z(J))*SINFHI(J)-(YTP-Y(J))*C0SPHI(J)1**2) 
DO 127 I = 1 , N1 

R ( I , J) = SORT ( (XT P-X ( I ) ) * *2 + (YT P-Y ( J) ) * *2 ♦ (ZTP-Z( J) ) **2) 
127 HL (I, J) = SORT ( (X(I )-XTP) ** 2 *■ HYZ(J)**2) 

VXC = 0. 0 
VYC = 0..0 
VZC = 0. 0 
C 

C CYCLE THROUGH VO'RTEX STRENGTHS, 

DO 150 K = 1, NN 
DO 150 L = 1,LL 
C 

C SELECT PARAMETERS FOR THIS PARTICULAR VORTEX STRENGTH, 

3 = L 
E = K + l 

IF (OPT1 ) GO TO 110 
A = L-l 

C = LL*2-L 

D = C-l 

IF (L-l) 150,129, 125 
125 IF (LL-L) 150,129 ,128 

110 IF (L-l) 150,113,111 

111 IF (LL-L) 150, 114, 112 

112 A = L-l 

C = MM-A 

D = MM-B 

GO TO 128 

113 A = MM 

C = MM 

0 » MM-1 

GO TO 128 

114 A = LL-1 


c 

617 

c 

618 

c 

619 

c 

620 

c 

621 

c 

622 

c 

623 

c 

624 

c 

625 

c 

626 

c 

627 

c 

628 

c 

629 

c 

630 

c 

631 

c 

632 

c 

633 

c 

634 

c 

635 

c 

636 

c 

637 

c 

638 

c 

639 

c 

640 

c 

641 

c 

642 

c 

643 

c 

644 

c 

645 

c 

646 

c 

647 

c 

648 

c 

649 

c 

650 

c 

651 

c 

652 

c 

653 

c 

654 

c 

655 

c 

656 

c 

657 

c 

658 

c 

659 

c 

660 

c 

661 

c 

662 

c 

663 

c 

664 

c 

665 

c 

666 

c 

667 

c 

668 

c 

669 

c 

670 

c 

671 

c 

672 
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C = LL+1 

C 

673 


0 = LL 

C 

674 


GO TO 128 

C 

675 

12 8 

RKA = R{ K, A) 

C 

676 


RKC = RC K »C)' 

c 

677 


REA = R( E . A) 

c 

678 


REC = R( E »C) 

c 

679 


HLKC = HL (K.C) 

c 

680 


HLEC = HL (E. C) 

c 

681 


HOA = HO ( A) 

c 

682 


HOC = HO (Cl 

c 

683 


YA = YCA) 

c 

684 


ZA = Z (A I 

c 

685 


ZC = Z(C> 

c 

686 


HYZA = H YZ (A) 

c 

687 


HYZC = HY Z (Cl 

c 

688 

129 

SINL = SINPHICU 

c 

689 


COSL = COSPHKL) 

c 

690 


RK8 = RCK.B) 

c 

691 


RKO = R(K,D)' 

c 

692 


RE 6 = R(E,8) 

c 

693 


REO = R( E . 0) 

c 

694 


HL KB = HL (K, 9) 

c 

695 


HLEB = HL (E. 91 

c 

696 


HOB = HO (B) 

c 

697 


HOO = HO ( 0) 

c 

698 


SIOEB = SIDE ( B ) 

c 

699 


OK = S(K) 

c 

700 


YB = Y (9 > 

c 

701 


ZB = Z(B ) 

c 

702 


ZO = Z(0) 

c 

703 


XK s X (K I 

c 

704 


XE = X(E) 

c 

705 


HY ZB = HYZ(B) 

c 

706 


HYZD * HYZ(D) 

c 

707 


P = 25.13274 

c 

708 

C 


c 

7 09 

C COMPUTE VELOCITY INDUCED BY VORTEX RECTANGLE OR RECTANGLES. TAKE ANY 

c 

710 

c special CASES INTO ACCOUNT. 

c 

711 


IF (L-l) 150,115,131 

c 

712 

131 

IF (LL-L) 15 0,115,132 

c 

713 

115 

IF (0PT1I GO TO 132 

c 

714 

13 0 

VXPS =0.0 

c 

715 


VYPS = 0.0 

c 

716 


VZPS =0.0 

c 

717 


IF (YTP.EQ.3.0) GO TO 230 

c 

718 


VXPS = l./(P'*SI0EB)*(HYZB*((RK0*RKB)*(SI0EB**2-( RKO-RKB) **2)7 ( ( 

c 

719 


1 HL KB* *2) *RK0*RK3) - (RED*RE3) * (SI OEB* * 2-( RE0-RE9 ) ** 2) / ( (HLE 3** 2) 

c 

720 


1*RE0*RE3'I ll*GAMA(K»H 

c 

721 

29 0 

IF (COSL. EO. 0.0) GO TO 66 

c 

722 


VYPS = (COSL/ (P’SIOEB) *(-( (RKO*R< B)* (SICE S**Z-(RKD-RK8)**2)/( ( 

c 

723 


2HLKB**2) *RK0*RK3) )*(XK-XTP) ♦ ( (RED+REB) * (SI0E9**2 -(RED-REB)** 

c 

724 


22) / ( ( HLE B**2) *RE0 *RE9) )* ( XE-XTP) ) * 1 ./ (P *DK) * ( ( ( RKB* RE 9) * ( DK** 

c 

725 


22- (RKB-RE B)**2)/( (HOB* *2 ) *RK8*RE3 ))*(Z9-ZTP)-((RK0*RED)*(DK**2- 

c 

726 


2 (RKO- REO) **2) / ( (HOO* *2) * RKO* REO) ) * ( ZO-ZTP ) I ) * GAM A (K ,L ) 

c 

727 


GO TO 67 

c 

728 
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65 VYPS - (l./(P*DK)*(( (RKB*REB)*(DK** C 

22- (RKB-RE B)**2)/( (HDB**2) *RKB*R£3 ) ) *( ZB-Z TP) - ( (R KD*REO) ♦ ( OK** 2- C 

2 (RKO-REO) **2)/((H00**2)*RK0*RE0) ) •( ZD-ZTP) 1 ) * GAM A (K,L ) C 

67 if cytp. EQ. a. a) go ro 231 c 

IF (ZTP.EQ.a.O) GO TO 201 C 

VZPS = (l./(»*OK) *(( (RK04-RE0) *(OK**2 -(RKO-REO)**2)/( (H0D**2) *RKD* C 
3R£0> -(RKB *■ REB)*(DK**2- <RK B-REB >**2 )/ ( ( HDB**2) *RKB*REB))* C 

3(Y8-YTP) ) )*GAMA(K,L) C 

201 VXC = VXC ♦ VXPS C 

VYC = VYC ♦ VYPS C 

VZC = VZC ♦ VZPS C 

GO TO 15 0 C 

13 2 VX = 0.0 C 

VY = 0.9 C 

vz = o.a c 

IF (YTP.EQ.0.0) GO TO 202 C 

VX s(i./(P*SIO£8>*((HYZ3*((RKA*RKB>*(SIDEB**2 -(RKA-RK3)**2>/ f ( C 

1HLK3**2) *RKA*RK-9) - (REA * RE3) * (SI DE B**2 - (REA-RE9)**2>/ ( (HLE3**2> C 

1*REA*RE9))) ♦ (HYZCM ( RKD *RKC ) *('S IDEB**2 - (RKO-RKC) ** 2) / ( ( C 

1HLKC**2) *RKC*RKO) - ( RED ♦ REC ) * (SI D£ 8**2 - (RED-REC>**2)/ ( (HLEC**2) C 
l*R£C*REO ) ))) >*GAHA(K,U C 

23 2 IF (COSL.EQ. 0.0) GO TO 58 C 

VY=(C0SL/<P*SI0E9)*(“((RKA*RKe>*(SID£3**2“(RKA-RKB)**2>/(( C 

2HLK3**2) *RKA*RK3) + (RKD+RKC ) * (SI DEB**2 - ( RKC-RKO ) **2 )/ ( (HLKC**2> C 

2*RKC*RK0> )*(XK-XTP) *■ ( ( RE A*RE8) * (SIOE 9** 2 -( REA-RE9) **2> / ( ( C 

2HL£B**2) *REA*RE3) ♦ (RE3*REO> * (SI DE8**2 - (REC-RED)**2)/ ( (HLEC**2) C 

2*REC* RED ) ) * ( XE-XT P) ) ♦ 1./ (P*OK> * ( ( (RK9 *R£B ) * ( OK**2 - ( C 

2RKE-REB) ** 21 / ( <HD0**2) *RK9*RE3)) * (ZB-ZT?) - { (RK9* RED) * (CK**2-(RKD- C 
2RE D)**2)/((H 00**2 ) *RKO*RE 0) )* (ZO- ZIP) ♦ ( (RKC*R'EC) * (DK**2- (RKC- REC) C 
2**2) / ( (HOC** 2) *RK C*R£C) ) * (ZC“ ZT<») » ( (RKA4-REA) * (DK**2-( RKA-REA) **2) C 
2/ ( (HO A** 2 ) *RKA*R£ A)) MZA-ZTP) ))*GAMA(K>L) C 

GO TO 69 C 

60 VY s ( l,/(F*OK)*(( (RK3+RE9)*( 0K**2 -( C 

2RKB-REB) **2) / ( (HDB**2) *RK B*RE B) ) * (ZB-ZT*) - ( (RKO*RED) * (OK**2-(RKO- C 
2REO)**2)/ ((H00**2)*RK0*RE0))* (ZO-ZTP) *((RKC*REC) ♦ (DK**2- ( RKC- REC) C 
2**2)/ ((HOC** 2) *RK C*REC) ) *(ZC-ZTP) - ( (RKA+REA)* (OK* *2- (RKA-REA) **2) C 
2/ ( (HOA** 2 ) *R!KA*RE A)) * { ZA- ZTP) ) )*G AMA (K» L) C 

69 IF (YTP.EO.O.O) GOTO 71 C 

IF (ZTP.EQ.0.0) GO TO 71 C 

IF (SINL.EQ. 0.000Q0) GO TO 70 C 

VZ * (SINL/(P*SIOE3)*(((RKA»RKB)* (SI0EB**2 -( P.KA-RKB) **2) / < ( C 

3HLKB**2) *RKA*RK3) - (RKD* RKD) * (SI DEB* *2 - (RKC-RKD)**2)/( (HLKC**2) C 
3*RKC*RK0) )*(XK-XTP) * ((REC*RED) * (SI0E3** 2 -(REC-REO) **2)/ ( ( C 

3HLEC**2) *RE3*REO) - (REA*RE9> * (SI DEB**2 - (REA-RE3)**2)/( (HLEB**2) C 
3*REA*RE3) )*(XE-XTP)) + 1. / (P*0K> * ( ( (RKAtREA)* (0K**2 - (P.KA C 

3-REA)**2)/((HOA**2)*RKA*REA) - (RKC+REC)* (OK**2 - <RKC-REC>**2)/ ( ( C 
3 HOC** 2)* RKC* REC) ) * (YA- YTP ) ♦ ( ( RKO*REO) * (0K**2 -(RKO-REO) **2)/ C 

3((H00**2) *RKO*REO) - (R<B*REB)*OK**2 - (RKB-RE3) * *2 )/ ( (HD 9**2) * C 
3RKB*REB) )*(YB-YTP)))*GA1A (K,L) C 

GO TO 71 C 

79 VZ = (l./(P*OK)*(((RKA«-REA)*(OK**2 - (RKA C 

3- RE A) **2) / ( (HOA** 2 )*RKA*REA) - (RKC+REO* (DK**2 - (RKC“R£C>**2)/ ( ( C 
3 HOC** 2)* RKC* REC) ) MYA-YTP) ♦ ( { RKD+REO) *(DK»*2 -( RK O-REO) ** 2) / C 

3 ( ( HDO**2 ) *RKO*REO ) - ( RKB *REB ) *( 3 K**2 - (RKB-RC3) **2)/ ( (HO B**2 ) * C 

3RKB*REB) ) * (VB-YTP ) ) ) *G AKA (K»L ) C 

71 VXC = VXC ♦ VX C 
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VYC = VYC «■ VY C 785 

V ZC = VZC ♦ VZ C 786 

150 CONTINUE C 787 

C C 788 

C COMPUTE VELOCITY INOUCEO BY MOOEL. C 789 

RM1 = RM (1) C 790 

RM2 s RM(2) C 791 

VXM = 0.0 C 792 

VYM = 0. 0 C 793 

VZM =0.0 C 794 

IE (HM1.LT. 1. E-10 ) GO TO 155 C 795 

VYM = GAMAM*2./P* (l.MXrP-XM>/RHl)/(HMl**2) C 796 

VZM = -VYMMYTP-YM) C 797 

VYM = VYMMZTP-ZNi) C 798 

155 I p (HM3.LT. 1. E-10) GO TO 160 C 799 

VXM = GAMAM*2./P* (1. «■ < XT P -XM) /RM2 )/ (HM3** 2) C 800 

VZM = VXMMYTP-YM) VZM C 801 

VYM = VXMMZM2-ZTP) *■ VYM C 802 

VXM = 0.0 C 803 

160 IF (HM2. LT.l.E-10) GO TO 165 C 804 

VXM = GAMAM* (RMlfrRM2) M5P AN** 2-(R M1-RM2) * *2) 7 (P* SPAN* RM1*RM2* C 805 

1 (HM2**2) ) C 806 

VYM = VYM + VXMMXM-XTP) C 807 

VXM = VXMMYTP-YM) C 808 

C C 809 

C COMPUTE TOTAL VELOCITY COMPONENTS. C 810 

155 VXT = VXC ♦ VXM C 811 

VYT = VYC *■ VYM C 812 

VZT = VZC ♦ VZM C 813 

RETURN C 814 

END C 815 
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ooooooooooooooooooooooo o o onooooooooonooooooooooooo o o o o n 


APPENDIX 0 


PROGRAM TO COMPUTE NON-LINEAR MIND TUNNEL WALL INTERFERENCE 
FACTORS FOR HIGHLY LOADED LIFTING SYSTEMS 


PROGRAM WINGT (IN PUT > OUT PUT, T APE5 = INPUT»T APE6-0UTPUT) 0 0 

0 1 

THIS PROGRAM IS WRITTEN IN FORTRAN I V FOR THE COO 6400 COMPUTER. 0 2 

APPROXIMATE STORAGE REQUIREMENT IS 52000 (OCTAL). D 3 

EXECUTION TIME IS APPROXIMATELY 230 SECONDS PER CASE WITH 29 TRAILING Q 4 

SEGMENTS, 7 ITERATIONS, ANO 100 SURREY POINTS. O 5 

D 6 

THE WINO TUN NEL CROSS-SECTION MUST HAVE A PLANE OF LATERAL SYMMETRY D 7 

AND MUST REMAIN CONSTANT OVER THE LENGTH OF THE TUNNEL O 8 

D 9 

INPUT DATA SEQUENCE O 10 

D 11 

I (ID D 12 

AN INTEGER PARAMETER WHICH DETERMINES THE Z COORDINATE OF TOP D 13 

ANO BOTTOM CENTER CONTROL POINTS. IF I.NE.l THESE CONTROL POINTS O 14 

WILL BE LOCATEO ON THE CENTERPLANE OF THE TUNNEL (I.E. Z=0.0>. D 15 

IF I.EQ. 1 THESE CONTROL POINTS WILL BE LOCATED AT Z(l)/2 O 16 

D 17 

MM, NN (212) D 18 

MM IS THE NUMBER OF COORDINATE PAIRS DEFINING THE COMPLETE CROSS- D 19 

SECTIONAL SHAPE OF THE TUNNEL. MM CANNOT EXCEED 20. D 20 

NN IS THE NUMBER OF VORTEX RECTANGLES MAKING UP THE LENGTH OF D 21 

THE TUNNEL. NN CANNOT EXCEEO 14. O 22 

D 23 

Y(I) , 7(1) (2F10.5) O 24 

Y ANO Z ARE THE COORDINATES, IN FEET, OF THE POINTS DEFINING THE O 25 

CROSS-SECTION SHAPE OF THE TUNNEL. MM CARDS ARE REQUIRED. D 26 


THE ORIGIN OF THE COORDINATE SYSTEM IS TAKEN ON THE TUNNEL CENTER O 27 
LINE WITH X POSITIVE OOWNSTREAM, Y POSITIVE UPWARD, AND 2 POSITIVE D 28 

TO THE RIGHT LOOKING DOWNSTREAM. THE FIRST CARD IN THE SEQUENCE IS D 29 

THE FIRST COORDINATE TO THE RIGHT (POSITI'VE Z> OF THE POSITIVE Y O 30 

AXIS, ANO SUBSEQUENT POINTS ARE TAKEN CLOCKWISE AROUND THE TUNNEL. O 31 


SEGMENT LENGTHS BETWEEN ADJACENT FOINTS SHOULD BE EQUAL, EXCEPT O 32 

THAT, IF CONVENIENT SPACING REQUIRES POINTS ON TOP ANO BOTTOM O 33 

CENTER LINE, THOSE POINrS ARE OMITTED AND THE FIRST DATA CARD D 34 

ABOVE, I, IS SET TO 1.0. O 35 

O 36 

OELTAX (F10. 5) 0 37 

LENGTH IN FEET OF THE VORTEX RECTANGLES IN THE STREAMHISE O 38 

DIRECTION. SHOULD BE EQUAL TO THE LENGTH OF SEGMENTS IN THE O 39 

CROSS-SECTION. O 40 

D 41 

BVOATA (F10. 5) D 42 

THE VORTEX SPAN OF THE WING IN F*EE AIR WHICH PRODUCED THE PUNCHED D 43 

CARD DATA TO 1 BE USED IN THIS PROGRAM.! 0 44 

D 45 

8V0TW (FI 0.5) 0 46 

THE RATIO OF VORTEX SPAN TO MAXIMUM TUNNEL WIDTH TO BE USEO IN 0 47 

THIS COMPUTATION. 0 48 

0 49 

YW(l) (F10.5) 0 50 
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C THE VERTICAL LOCATION OF THE MODEL BOUND VORTEX IN THE TUNNEL. 

C 

C OELTAX (F10.5) 

C THE VORTEX SEGMENT LENGTH TO BE USED IN CONSTRUCTING THE 

C TRAILING VORTEX PAIR IN THE TUNNEL. NEED NOT BE THE SAME 

C AS THAT USED IN THE FREE AIR PROGRAM, USUALLY SPAN/10 

C 

C ZMAX, YMIN C2F10.5) 

C MAXIMUM Z COORDINATE AND MINIMUM Y COOROINATE TO BE ALLOWED 
C IN SURVEY OF WALL INTERFERENCE VALUES. THESE PARAMETERS WILL 

C BE USED TO DETERMINE IF A SURVEY POINT LIES TOO NEAR THE TUNNEL 

C FLOOR OR SIOE WALLS FOR ACCURATE INTERFERENCE COMPUTATION. 

C 

C SPAN, SPEED, GAMAM, ASPECT, FAL, VXWC , VYWC, ALFA'. (4E20. 10) 

C THESE THREE CAROS OEFINE THE MODEL TO BE USED IN THIS COMPUTATION, 
C AND ARE PART OF THE OECK PUNCHED BY THE WING-IN-FREE- AIR PROGRAM. 

C SPAN IS WING VORTEX SPAN, FEET. 

C SPEED IS REMOTE WIND SPEED IN THE TUNNEL, FEET/SECOND 

C GAMAM IS MODEL WING CIRCULATION, FEET SOU AREO/SECOND. IF GAMAM IS 

C LESS THAN OR EQUAL TO ZERO, THE ZERO LIFT CASE IS PERFORMED. 

C ASPECT IS THE ASPECT RATIO OF THE WING. 

C FAL AND F AO ARE THE LIFT AND CRAG OF THE WING IN FREE AIR, POUNDS. 

C VXWC AND VYWC ARE VELOCITY COMPONENTS AT THE CENTER OF THE SOUND 

C VORTEX IN FREE AIR. 

C ALFA IS THE WING ANGLE DF ATTACK IN FREE AIR. 

C 

C XFA, YFA, ZF A , VXTOT, VYTOT, VZTOT. (4E2Q.10) 

C THESE ARE THE COORDINATES ANO VELOCITIES SURVEYED BY THE 
C WING-IN- FREE-AIR PROGRAM AND PUNCHED IN A. CARD DECK. THE 
C COORDINATES ARE REFERENCED TO THE WING. 

C NOTE THAT ZF A AND YFA ARE ALSO PROGRAM BRANCHING PARAMETERS. 

C IF (ZFA. EO. 10000. ) THE PROGRAM TRANSFERS TO NEW MODEL OATA. 

C THEN IF (Y FA .EO. 10003. ) THE PRESENT MODEL WAKE COORDINATES 

C ARE USED FOR THE FIRST ITERATION OF THE NEW WING. THIS REDUCES 

C THE NUMBER 0!F ITERATIONS. 

C 

1 FORMAT (212) 

2 FORMAT (2F10.5) 

3 FORMAT CF10. 5) 

4 FORMAT (4F13.5) 

5 FORMAT (1F10.5I! 

6 FORMAT (12) 

7 FORMAT ( 2F10.5) 

8 FORMAT (ID 

q FORMAT (I3,7F10.5) 

11 FORMAT (4E20.1U) 

12 FORMAT (3F10.5) 

13 FORMAT (5F1J.5) 

30 FORMAT ( 1H1, 13X,19HTUNNEL COOROI NATES,/, ✓, 14X, 1 HY, 13X , 1H Z, 

1 (/,10X,F10.5,4X,F10.5)) 

31 FORMAT (/ ,/, 15X,1 OHX STATIONS, (/ , 4X, 5F6. 2 ) ) 

32 FORMAT ( 1HQ, 22HCR0SS SECTIONAL AREA = ,F10.4> 

5300 FORMAT (1H0.14HTAIL LENGTH = ,F7. 2,5X , 14HTAIL HEIGHT = ,F6.2,5X, 

1 19HSPANW I SE STATION = ,-6.2) 

5010 FORMAT ( 1H ,13H(F.A.) VX = , F9. 3 ,6X, 5HVY = , F9. 3 ,6X, 5HVZ = , 
1F9.3,6X, 8H ALPHA = ,F7. 4, 6 X,7HBETA = ,F7.4> 
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D 

83 

0 

84 

D 

85 

0 

86 
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50 20 FORMAT (1H , 13H(TUN.) VX = , F9. 3 ,6X, 5HVY = , F9. 3 ,6X, 5HVZ = , 0 

1F9.3,6X,8HALPHA = ,F7.4,6 X,7HeETA = ,F7.4) 0 

50 30 FORMAT CiH , 13H(C0R.) VX = , F9. 3 ,6X,5HVY = ,F9. 3,6X, 5HVZ = , 0 

1F9.3,6X, 6HALPHA = ,F7. 4, 6 X, 7HBET A = ,'F7.4> 0 

5040 FORMAT <1H , 34HCORRECTION FACTORS DEL (ALPHA) = ,F8. 3, 10X, 12HDEL ( D 
13ETA) = ,F8. 3,10X,5HOQ = ,F8.4) 0 

5100 FORMAT ( 1H1, 12HWING SPAM = , F6.2, 1CX, 8HGAMAM = , F7. 2, 10X , 15HASPECT D 
1 RATIO = ,F5. 2,10X,18HREMOTE VELOCITY = >F8.2) D 

5110 FORMAT (1H ,23H(F.A. CENTER) LIFT = , F8 . 3,10X, 7H0RA6 = ,F7.4, 0 

110X,5HVX = ,F8.3,1CX,5HVY = F8.3) 0 

5120 FORMAT <1H ,23H(TUN. CENTER) LIFT = , F8 .3, 1QX, 7 HORAG = ,F7.4, D 

110X,5HVX = , F8.3, 10X,5HVY = ,F8.3) D 

5130 FORMAT CIH , 23HCC0R. CENTER) LIFT - ,F8 . 3, 10X , 7H0RA G = ,F7.4» D 

110X,5HVX = , F8.3, 10X,5HVY = ,F8.3) 0 

5140 FORMAT (1H » 34HCORRECTION FACTORS' DEL(ALPHA) = ,F8. 3,10X, 0 

15H0Q = ,F8.4) 0 

DIMENSION XC15),YC20) , Z( 2 0) ,SINPH I< 20 ) ,COSPHI (20 ) ,XCPT(14>, D 

lYCPT(lO) , ZCP’T (10) ,R(15,20) ,SI0E(2C> ,HL(15,20) ,HD ( 20) , S (14) , ZM ( 2) , 0 

1HM(3) ,HYZ(20) ,RM(2),GL(10 ) 0 

DIMENSION CC (100, 109) ,GAM A(14,1Q> ,GAMAK(1C0,1) 0 

DIMENSION XW(4C> , YW(4G> , Z W(40) ,RW(2,2) ,DSM(39) ,V 0AR(2 > D 

LOGICAL STWK,0PT1 0 

REAL LIFT 0 

RHO a .002378 D 

YFA=0.0 0 

14 CONTINUE 0 

C 0 

C READ DATA DESCRIBING TUNNEL FROM PUNCHEO CARDS. D 

REAO (5,3) I D 

OPT 1 * I.EQ.l 0 

34 REAO (5,1) MM,NN 0 

REAO (5,7) (Y(I),Z(I), 1=1, MM) D 

REAO (5, 3) OELTAX 0 

C 0 

G TEST DIMENSIONS 0 

' IF ( ( MM, GT .20). OR. (NN.GT. 14) ) GO TO 906 0 

C D 

C TEST SCALING OF TUNNEL, IF NECESSARY CHANGE SCALE SO THAT THE WING 0 

C SPAN OF MODEL IN TUNNEL CORRESPONDS TO THAT OF MODEL IN FREE AIR. D 

XCI = Z(l) D 

C READ SCALING DATA FROM PUNCHED CAROS. 0 

READ (5,3) 3VDATA D 

READ (5,3) 3V0TW D 

DO 35 I = 2, MM 0 

IF (Z (I) . GT. XCI) XCI = Z(I) 0 

35 CONTINUE D 

YCJ = BV0ATA/8V0TW/2. D 

XCI a. YCJ/XCI 0 

C IF THE SCALING FACTOR IS UNITY DO NOT CHANGE TUNNEL SIZE. 0 

IF (XCI.EO.l.) GO TO 37 D 

DO 36 1=1, MM D 

V(I> a Y (I)*XCI 0 

Z(I) = Z(I)*XCI D 

35 CONTINUE D 

OELTAX = OELTAX*XCI D 

37 CONTINUE 0 
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c 

C COMPUTE THE CROSS SECTIONAL AREA OF THE TUNNEL. 

AREA =0.0 
00 38 I = 2, MM 

AREA = AREA * ABS (Y(I) -Y ( 1-1) )*ABS (Z tl» *Z (I-l>> 

33 CONTINUE 

AREA = AREA/2. 

C 

C 

C NOW COMPUTE THE TUNNEL PARAMETERS TO BE USED IN THE COMPUTATION. 

LL = MM/2 ♦ 1 
NL = NN * LL 

IF (NL.GT.100) GO TO 906 

NM = NN* MM 

Nl = NN + 1 

XCI = 0.0 

DO 20 1=1, NN 

X(I) = XCI 

20 XCI = XCI ♦ DELTA X 

X ( Nl) = 1000.0 «- X (NN) 

DO 21 I = 1, NN 

21 S(I)=X(I*i) - X<I) 

00 23 1=2, MM 

22 SIOE(I) = SORT ( (Y (I) - Y(I-1))**2 * (7(1) - Z(I-1>>**2) 

SINPHKI) = ((Y(I)-Y(I-l) )/(SIDE(I>)) 

23 COSPHICI) = ((Z(I)-Z(I-l) )/(SIOE(I>>> 

SIOE ( 1) = SORT ( (Y (1) - Y ( MM) ) **2 ♦ (Z(i) - Z(MM))**2)' 

SINPHICl) = ((Y(1)-Y(MMI)/(SIQE(1>)) 

COSPHICI) = (CZ(l)-Z(MM))/(SIOE(l))) 

DO 24 I = 2, LL 

YCPT(I) = (Y (I)*Y (I-l))/(2.) 

24 ZCPT(I) = (Z(I)*Z(I-l)>/(2.) 

ZCPTC1I = CZ (ll + Z (MN) )/(2. ) 

YCPTCU = (Y (1)*Y <MM> )/(2.) 

IF (.NOT. OPT 1) GO TO 91 
ZCPT(l) = Z(l)/2. 

ZCPTILLI = Z(LL-l)/2. 

91 MMM=NN-1 

00 25 I = 1, MMM 

25 XCPT(I) = (X(I+1) ♦ X(I))/(2.) 

XCPT(NN) = X (NN) *■ DELTAX/2, 0 

C ALL TUNNEL PARAMETERS HAVE 3EEN COMPUTED.! 

C 

C 

C GENERATE THE MATRIX OF COEFFICIENTS. 

CALL MATRIX (MM, NN , LL , N 1,X,IY, Z, SINP HI ,C0SPHI, SI OE , S , XCPT, 

1 YCPT , ZCPT ,CC) 

C 

C INVERT THE MATRIX OF COEFFICIENTS. 

CALL INV R (CC,NL,1 00, 100) 

C 

C WRITE A DESCRIPTION OF TUNNEL. 

WRITE (5,30) (Y(I),Zm, 1 = 1, MM) 

WRITE (6,31) (X(I), 1=1, Nl) 

WRITE (5,32) AREA 
C 
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C READ MODEL INFORMATION FROM PUNCHED CAROS. 

READ (5,3) YW(i) 

READ (5,3) DELTAX 
READ (5,7) ZMAX, Y MIN 

15 CONTI NIE 

REAO (5,11) SPAN, SPEED , 5 A MAM , ASP'S CT » FAL»F AD, VXWC » VYWC » ALFA 
IF (EOF, 5) 907,16 

16 CONTINUE 

IF (YFA. EQ. 10000. ) GO T3 40 
C 

C NOW GENERATE MOOEL PARAMETERS. 

IF (GAMAM. GT. 0.0) NW=3Q 
I = NN/2 
XW(1) = X(I) 

XW (2) = XW(1) 

YW (2) = YU (1 ) 

ZW ( 1) = 0.0 
ZW (2) = SPAN/2.' 

ZH (3) = ZW(2) 

STWK= (GAMAM.LE.O. 0) 

IF (STWK ) GO TO 18 
NW1 = NW + 1 

CHORD = SPAN/ (ASPECT*. 785 398163** 2) 

ALFAA = A SIN (GAMA M/( 3. 141 5927 *CHO RD* SPEED) ) 

XCI = C. 75*CHORD*SQRT (1.- (.78539816**2)) 

XW ( 3) = X W (2 ) *■ XCI*COS( ALFAA) 

YW (3) = Y W (2 ) - XCI*SIN(ALFAA) 

XCI = OELTAX ♦ XW (3) 

YCJ = YW ( 3 ) 

7CJ = ZW(3) 

DO 90 N = 4, NW 
ZW(N> * ZCJ 
YW (N) = YCJ 
XW (N) = XCI 
XCI = XCI ♦ OELTAX 
93 CONTINUE 

XW (NW 1) = XW (NX) + 1000.0 
YW(NWl) = YCJ 
ZW(NWl) = ZCJ 
GO TO 19 
C 

C IF THE STRAIGHT WAKE (ZERO LIFT COEFFICIENT) SOLUTION IS REQUIRED 
C SET UP A HORSESHOE VORTEX MODEL. SET SPEED TO 1009., GAMAM TO 1.0. 

IS XW (3) = XW (2 ) * 1000. 

YW ( 3) = YW(2) 

S°EEO = 1000. 

GAMAM = 1.0 
C 

C COMPUTE THE LIFT AND INDUCED DRAG OF THE WING IN FREE AIR. 

FAL = RHO*SPEED*SPAN*G AM AM 
FAD = RHO* (G AMAM* *2) / 3. 14 159 
NW = 2 
NW1 = NW ♦ 1 
19 OO 81 I = 1, NW 

J = 1*1 

81 DSM(I) = SQRT((XH(I)-XW(J))**2*(YW(I)-YW{ J))**2* (ZW(I)-ZW(J))**2) 
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Q = . 5*RH0*SP£E0**2 
C 

C BEGIN ITERATIVE PROCEOURE. NUMIT IS THE NUMBER OF ITERATIONS TO BE 
C USEO. IF THIS CASE REPRESENTS A SMALL CHANGE FRO'M A PREVIOUS 
C EQUILIBRIUM STATE, REOUCE NJMIT. 

NUMIT = GAMAM/19. «■ 2. 

40 CONTINUE 

IF (YFA.'EQ.iOOOO. ) NUMIT = GAMAM7 30. 2. 

900 DO 901 NUMBER = 1, NUMIT 
C 

C COMPUTE THE RIGHT HANO SIOE OF THE MATRIX EQUATION. 

CALL RHS CXCPT.YCPT ,ZCPT,XW ,YW,iZW,DSM,GAMAM ,SPAN,SPEEO, 

1GAMAK,NN,LL»NW,SI NPHI , COS PHI ) 

c 

C COMPUTE THE VORTEX STRENGTHS. 

1001 00 101 1 = 1, NL 

J = (I-1)/LL ♦ 1 
K = (1-J)*LL f I 
XCI = 0. 0 
00 100 L = 1 , NL 

100 XCI = XCI4-CC (I, LI *GAMAK(L ,1) 

GAMA ( J,K ) = XCI 
131 CONTINUE 
C 

C IF THIS IS WITHIN THREE ITERATIONS OF THE LAST WRITE COMPUTED VORTEX 
C STRENGTHS. 

IF ( ( NUN I T-N UMBER) *GE. 3) GO TO 110 

3999 FORMAT ( 19H1ITERATI0N NUMBER ,12) 

WRITE (6,3999) NUMBER 

50 0 FORMAT ( 3 OHO CALCULATE!) VORTEX STRENGTHS ) 

WRITE (6,503) 

00 502 J = 1 , NN 

WRITE (6,501) (GAMA( J, K) , K=1,LU 
501 FORMAT ( / , 11F11. 5 ) 

50 2 CONTINUE 

4000 FORMAT ( 27H1 RESULT ANT V3RTEX STRENGTHS ) 

40 02 FORMAT (13H3RING NUM3ER , 12,/, 11 F11.-5) 

4004 FORMAT ( 15H0 SECTI ON NUMBER ,13 ,/, llFll.5) 

4010 WRITE (6,4000) 

4315 00 4140 L=i, N1 

40 20 M=L-1 

4325 00 4375 I=1,LL 

40 30 IF (L-2) 4C50, 4060, 4040 

4040 IF (L-Nl) 4360, 4070, 4140 

4050 GL (I) = GAMAi(L,I) 

43 55 GO TO 40 75 

4360 GL ( I) = GAMA (L,I) - GAM1(M,I) 

43 65 GO TO 40 75 . , 

4070 GL (I) = -GAMA (M,I ) 

40 75 CONTINUE 

4080 WRITE (6,4002) L, (GL (I ) , 1=1, LL) 

4100 IF (L-Nl) 4110, 4140, 4140 
4110 00 4125 1 = 2, LL 

4115 J=I-1 

4120 GL (J) = GAMA (L, J) - GAMA ( L,I) 

4125 CONTINUE 
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*♦130 HMM = LL - 1 0 

4135 WRITE 16,4004) L, (GL( J) , J =1, MMH) 0 

4140 CONTINUE 0 

C 0 

C 0 

C PERFORM WAKE ITERATION PROCESS. 0 

110 CONTINUE 0 

CALL WKIT ( XW, YW, ZW , X, Y ,Z, S INPHI ,CO SPHI ,SI DE,S, GAMA ,OSM , 0 

1GAMAM, SPEED, SPAN, NW,NN, MM , Ni , LL,N W1 ,RHO,Q ,FAL, FAD, CHORD, LIFT, DRAG, D 
1STWK, VXTC,VYTC,ALPHAO,ALPHAI, ALFA'A,VXMC,VYMC) 0 

C D 

C IF THIS IS THE LINEAR CASE (ZERO LIFT COEFFICIENT) GO DIRECTLY TO D 

C WALL CORRECTION SURVEY, DO NOT PERFORM ANY ITERATIONS. D 

IF (STWK) GO TO 810 D 

901 CONTI NIE 0 

GO TO 811 D 

C D 

C COMPUTE VXWC ANO VYWC FOR S 3 ECIAL CASE OF ZERO LIFT COEFFICIENT. 0 

810 VXWC = VXMC ♦ SPEED D 

VYWC=VYNC 0 

C D 

C WRITE A DESCRIPTION OF MODEL ANO TUNNEL OPERATING CONDITIONS. 0 

811 WRITE (5’, 4240) GAMAM 0 

4240 FORMAT € 1H0, 19HM0 DEL CIRCULATION =,F10.5> 0 

WRITE (6,4195) SPAN D 

4195 FORMAT ( 14H0 VORTE X SPAN = ,F10.5> 0 

WRITE (6,4200) Q D 

4200 FORMAT ( 11H0TUNNEL Q = ,F10.5) D 

WRITE (6,4185) SPEED D 

4185 FORMAT ( 26H3TUNNEL NOMINAL VELOCITY s ,F10.5) D 

WRITE (6,5100) SPAN, GAMAM, ASPECT, SPEED 0 

C D 

C WRITE FREE AIR RESULTS, 0 

WRITE (6,5110) FAL, FAD, VXWC, VYWC 0 

C D 

C WRITE TUNNEL RESULTS. D 

WRITE (5,5120) LI FT, DRAG, VXTC ,VYT C 0 

FAL=FAL-LIFT D 

FAO=F AO- DRAG ' D 

OA= VXWC- VXTC D 

OB=VYWC-VYTC D 

C 0 

C WRITE CHANGES OUE TO TUNNEL.* D 

WRITE (6,5130) FAL, FAD, OA ,D9 D 

OA=(ATAN(-VYWC/VXWC>-AT4N(-VYTC/VXTC))*AREA*Q/LIFT 0 

DQ= (VXWC* *2* VYWC* *2- VXTC* *2- VYTC* *2)7 (SPE ED** 2) D 

C D 

C WRITE ANGLE OF ATTACK CORRECTION FACTOR ANO DYNAMIC PRESSURE RATIO. D 
WRITE (6,5140) 0A,00 D 

801 CONTINUE D 

REAO (5,11) XFA,Y FA,ZFA,VXTOT ,VYT OT ,VZTOT D 

IF (EOF, 5 ) 907,802 D 

80 2 CONTINUE 0 

IF (ZFA. EQ.1000Q. ) GO TO 15 0 

IF (ZFA, GT.ZMAX) GO TO 801 0 

OA=ATAN( YFA/XFA) D 
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DB=ALFA*0A 

TL=SQRT( YFA**2*XFA**2) 

TH=TL*SIN(03) 

TL=TL*COS (D3) 

XCI=XW(1) ♦TL*COS( ALFAA)*TH*SIN(ALFAA) 

YCJsYW(l)-TL*SIN( ALFA A) ♦ TH*C0S(ALFAA) 

ZCJaZFA 

IF (YCJ.LT.YMIN) GO TO 801 

CALL XYZVEL (XCI, VC J, ZC J, XW ,YW, ZW,X , Y ,Z ,SI NPHI, COSPHI ,SIDE , S 

1,GAMA,0SM,GAMAM,SPE£9,SPAN,NW,NN, MM,N1,LL ,VXT,VYT,VZT ,VXR,VYR, 
1VZR,VXM,VYM,VZM) 

IF (• NOT • STWK) GO TO 853 

VXTOT=VXM*SPE£D 

VYT0T=VYM 

VZT0T=VZM 

850 WRITE (6 , 5C0 0) TL,TH,ZCJ 
AL PHA=AT AN (** V YTOT/ VXTOT) 

3ETA=ATAN(VZTOT/ VXTOT) 

WRITE (6,5010) VXTOT ,VYTOT ,VZTOT, ALPHA, BETA 
OA=ATAN(-VYT/VXT) 

OB=ATAN( VZT/VXT) 

WRITE (6,5020) VXT, VYT ,VZT, DA ,09 

DQ=VXT0T**2*VYT0T**2*VZT0T**2 

VXTOT=VXTOT-VXT 

VYTOT=VVTOT-VYT 

VZT0T=V7T0T-VZT 

OA=ALPHA-DA 

OB=BET A-OB 

WRITE (5,5030) VXTOT , VYT 0 T , VZTOT, DA ,09 

OA=OA*AREA*Q/LIFT 

0B=D3-AREA*0/LIPT 

OQ=(OQ-( VXT**2*VYT**2*VZT**2> )/(SPEEO**2) 

WRITE (6,5040) OA,OB,OQ 
GO TO 801 

906 CONTINUE 
WRITE (6,902) 

932 FORMAT (52H0 DIMENSIONED STORAGE EXCEEDED - EXECUTION TERMINATED ) 

907 CONTINUE 
STOP 
END 
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SUBROUTINE MATRIX (MM , NM , LL, Ni, X,;Y, Z, SINP HI ,COSPHI , SI OE,S ,XCPT , 0 

1YCPT,ZC?T,CC:> 0 

C 0 

C THIS IS A SUBROUTINE TO COMMUTE THE MATRIX OF COEFFICIENTS. 0 

C 0 

DIMENSION X(15)»Y (20) ,Z(20),SIN o HI(20),C0SPHI(20)»XCPT(i4)* D 

lYCPT(iQ) , ZC°',T ( 13) , R( 15 *20 ) ,SIDE(20) ,HL ( 15 ,20) , HO (20 ) , S (14 > , ZM C 2) , 0 

1HM (3) *HYZ(2Q) *RM(2)*GL(10) D 

DIMENSION CC(109, 100) D 

INTEGER A * 8* C * D*E 0 

M = 0 D 

C D 

C CYCLE THROUGH CONTROL POINTS (I.E. ROWS OF COEFFICIENT MATRIX). D 

DO 50 I * 1,NN 0 

00 49 J = 1,LL D 

M = M * 1 0 

C D 

C D 

C RECALL PARAMETERS FOR THIS CONTROL POINT,t GENERATE PARAMETERS FOR D 

C VORTEX NET WITH RESPECT TO THIS CONTROL POINT. D 

P = 25.13274 D 

SI N J * SINPHI(J) D 

COSJ = COSPHI(J) D 

XCI = XCPT(I) D 

YCJ = YCPT(J) D 

ZCJ = ZCPT(J) 0 

37 DO 26 JJ = 1,NM D 

HD(JJ) = SORT ( (YC J-Y(JJ) ) **2 + (ZCJ-Z ( JJ) ) **2) D 

HYZt JJ) = SQRrt ( (ZC J-ZCJJ) ) *SINPHI( JJ) - (Y C J-Y ( JJ ) ) *COSPHI ( J J) ) **2) D 
00 26 11*1, N1 0 

R(II, JJ) =SGRT< (XCI-X(II) ) **2f (YCJ-Y (J J) ) **2*( ZCJ- ? ( JJ) ) **2> D 

25 HL<II,JJ)*SQRT(CX(II)-XCI)**2 ♦ MYZ(JJ)**2) 0 

N = 0 D 

C D 

C CYCLE THROUGH VORTEX RECTANGLES (I.E. COMPUTE ELEMENTS IN THIS ROW D 

C OF THE COEFFICIENT MATRIX). D 

00 48 K= 1 ,NN D 

DO 47 L=1,LL 0 

N = N * 1 0 

C D 

C RECALL VARIABLES FOR THIS PARTICULAR VORTEX RECTANGLE PAIR. 0 

A a (L-l ) D 

B = L 0 

C * 2*LL-L D 

0 = C-l 0 

E = K*1 0 

C 0 

C IF THIS IS A SINGLE VORTEX ON THE TOP OR BOTTOM OF THE TUNNEL, NOT 0 

C ALL PARAMETERS ARE NEEDED. D 

IF (L-l) 50, 29,27 D 

27 IF (LL-L) 50, 29, 28 D 

28 RKA a R( K, A)’ D 

RKC = RCK,C) D 

REA a RCE,A) 0 

REC = R(E,C) D 

HLKC = HL (K, C) 0 
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O O O o O 


HLEC = HLCE.C) 

HOA = HO ( A) 

HOC = HO tC) 

YA = Y (A ) 

ZA = Z (A ) 

ZC = Z(C) 

HY ZA = HYZ(A') 

HYZC = HYZ(C> 

29 SINL = SINPHICL) 

COSL = COSPHI(L) 

RKB = RC K , B) 

RKO = RCK.D) 

RE8 = RCE,B)< 

REO = RCE ,0) 

HLKB = HL CK.B) 

HLEB = HL(E,8) 

HOB = HOC 8) 

HOO = HO (0) 

SI0E8 = SI0EC8I 
0K=SCK) 

YB = YCB) 

ZB = ZCB) 

ZD = ZCO) 

XK = XCK» 

XE = XCEI 
HY ZB = HYZC9) 

HYZD = HYZCOI 

IP (COSJ.EO.O. 000 00) GO TO 35 


COMPUTE THE Y, Z VELOCITY COMPONENTS INDUCED BY VORTEX RECTANGLE 
OR RECTANGLE PAIR. 

USE EQUATIONS APPLYING TO VARIOUS SPECIAL CASES. 

IF(L-l) 5C, 33,31 

31 IF CLL-L ) 50.33,32 

32 IF CCOSL.EQ. 0.G00Q0) GO TO 62 

VY= (COSL/ (P*SI0£9)*C-( (RK A*RK B)* C SIOEB**2- ( RK A-RKB) *• 2)/ C t 
2HLKB**2) *RKA*RK9) ♦ CRKO f RKC) * (SI OE B**2 - (RKC-RK D )**2 > / C ( HL KC**2 > 
2*RKC*RKO) )*CXK-XCI) ♦ (( R EA*RE8) * (SIDED** 2 -( PE A-RE B) **2) /( < 
2HLEB**2) *REA*RE3> ♦ (REC* REO) * (SI OEB**2 - (REC-RED>**2)/( CHLEC**2> 
2*REC*RE0) )*(!XE-XCI)> * 1. / (P*0< )* ( ( (RKB «-RE3> *( OK**2 -( 

2RKB-REB) **2)/ ( (HOB** 2) *RKB*REB)) * (ZB-ZC J> - ( (RKO* REO) * (DK**2-( RKD- 
2RED) **2) / ( CHOD**2 ) *RKD*RE O) ) * (ZO- ZC J> * ( (R KC*REC) * (OK** 2- ( RKC-REC ) 
2**2) / C (H0C**2)*RKC*REC>) * (ZC-ZCJ) -{ (RKA*REA) * (OK* *2- ( RKA-RE A) **2) 
2/ ( (H0A**2)*RKA*REA))*(ZA- ZCJ) )) 

GO TO 36 

62 V Y = (l,/(P*OK)*(( (RKB+REB)*(OK**2 -( 

2RKe-REe>**2) / ( (HO 8**2> *RK B*RE 9) ) * (ZB- ZCJ) - ( (RKO* REO) * (0K**2-(RKD- 
2RE0)**2) / ( (H00**2)*RK0*RE0)) * (ZO- ZCJ) t ( (RKC*REC) * (OK* * 2 - (RKC-REC) 
2** 2)/ C CHDC** 2) *RKC*REC) ) ♦ (ZC-ZCJ) -CCRKAfREA)* (OK**2-( RKA-RE A) **2> 
2/ ( CHO A** 2)*RKA*REA))*(ZA- ZCJ) )) 

GO TO 36 

33 IF CCOSL.EQ. 0.00000) GO TO 63 

VY = CCOSL/ (P*SIQ£B)*(-C (RKO*RK3) *(SI0E9**2-CRKD-RKB) **2)/C C 
2HLKB**2) *RKO*RKB) ) *CXK-XGI) ♦ (CREO*REB) * (SI0EB**2 -(REO-REB)** 
22) / ( C HLE B**2)*RE0*RE9) ) *C XE-XCI) ) ♦ l./(P*DK) *( (( RKB* RE 9) *(OK** 
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22-(RKB-R£8>**2>/( (H33**2 ) *RKB*RE9) ) M Z8-ZCJ)- ( (R KO«-RED) MDK**2- 
2(RKD-R£0I **Z) / i (H 00**2) *RKO* RED) ) MZ0-7CJ))) 

GO TO 36 

63 VY = (1./. (P*OK) *((l(RKB*REB> * (0<** 

22-(RK8-R£B>**2)/( (HDB**2> *RK B*RE8 > > * t ZB-ZCJ) - ( (RKDfRED) MDK** 2- 
2 (RKO-RED) **2)/ ( (H 00**2) *RKD*RED) ) MZO-ZCJ))) 

GO TO 36i 

33 VY = 0.30000 

36 IE (SINJ.EQ. 0.000 00) GO TO 42 

IF(L-i) 50,40,38 
39 IF(LL-L) 50,40,39 

39 IF (SINL.EO.O. 00000) GO TO 64 

VZ = (SINL/(P*SIOEB> *( ((RKA+RKB) * <SI0E8**2 -(RKA-RKB) **2>/ ( ( 
3HLKB**2) *RKA*RKB) - (RKC+RKO) MSI DEB**2 - CRKC-RK D)**2)/ ( ( HLKC** 2) 
3*RKC*RK0) )*('XK-XCI) ♦ ( ( REC*REO) * (SIOEB** 2 MPEC-RED) **2) / ( ( 
3HLEC**2) *REC*REO) - (REA*REB)*(SI0E8**2 - (REA-RE B) **2) / ( (HLEB**2> 
3*REA*RE3) ) M XE-XCI)) + 1. / (P*DK) * < ( (RKA+REA) * (OK **2 - (RKA 
3-REA)**2)/((H0A** 2)*RKA*REA) - (RKC+REC) * (0K**2 - (RKC-REC) **2) / ( ( 
3HDC**2)*RKC*REC)> MYA-Y3J) ♦ (( RKD+REO) MDK**2 -(RKD-REO) **2)/ 

3((H00**2)*RKD*RE0> - (R<B + REB)M0K**2 - (RK8-RE8) **2)/ ( (H0B**2) * 
3RKB*RE9) )MYB-YCJ)>> 

GO TO 43 

64 VZ = (l./(P*OK)* (((RKA«-REA)*(OK**2 -(RKA 

3-REA)**2)/(T'HDA**2)*RKA*REA> - (R KC+REC) • (OK**2 - (RKC-REC)**2>/ ( ( 
3HOC**2)*RKC*REC)> MYA-YU) ♦ ( ( RKD+REO) MDK**2 -(RKD-RE0>**2>/ 

3 ( ( HD0**2 ) *RK D*RE9 ) - (RK3*RE9)MDK**2 - (R KB-REB) **2>M (HD8**2> * 
3RKB*REB) 1 MYB-YCJ) >) 

GO TO 43 

43 VZ = (W(P*OK)M ((RK0ME0>*(0K**2 -(RKD-PED)**2) /( (HDD**2) *RKO* 

3 RED) -(RK8 ♦ REB) * OK**2- (RK9-RE3) **2>/ (<!HDB**2) *RKB* RE9) > * 
3(YB-YCJ))I 
GO TO 43 

42 VZ = 0.0 0 000 
C 

c 

C THE VELOCITY COMPONENTS HAVE BEEN COMPUTED, STORE' THE NORMAL VELOCITY 
C AT THIS CONTROL POINT IN CC ARRAY ELEMENT M,N. 

43 CC (M, N) = VY’COS J - VZ*SINJ 
C 

47 CONTINUE 

43 CONTINUE 

49 CONTINUE 

50 CONTINUE 
C 

C THE MATRIX HAS BEEN GENERATED, RETURN TO CALLING PROGRAM. 

C 

RETURN 

END 
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ooo ooo ooo o o o oo ooooooooo 


SUBROUTINE INVRCA ,N,ISIZE ,JSIZE) 0 587 

0 588 



SUBROUTINE TO COMPUTE THE INVERSE OF A MATRIX OF SIZE LESS THAN 

O 

589 


OR EQUAL TO 100 

O 

590 



D 

591 

THE 

MATRIX A IS REPLACEO BY ITS INVERSE. 

O 

592 

THE 

MATRIX IS ASSUMED TO CONTAIN N ROWS AND COLUMNS. 

O 

593 

ISIZE AND JSIZE ARE THE DIMENSIONS OF A. 

O 

594 

NOTE THAT THIS SU9R0U TINE DOES NOT TEST THE SINGULARITY OF A. 

O 

595 



D 

596 


DIMENSION IPIVOTt 1G0) , A<ISIZE,USI ZE > ,-INDE X (10:0 ,2 ) ,PIVOT{100> 

D 

597 


EQUIVALENCE (IRON, JROW), ( ICOLUM, J COLUM) , ( AMAX ,T, SWAP) 

O 

598 



D 

599 



O 

600 

15 

DO 20 J=1,N 

O 

601 

20 

IPI VOT (J ) =0 

D 

602 

30 

OO 550 1=1, N 

O 

603 



D 

604 


SEARCH FOR PIVOT ELEMENT 

O 

605 



D 

606 

40 

AMAX= 0.0 

D 

607 

45 

DO 105 J=1,N 

D 

608 

50 

IF CIPIVOT(J)-l) 60, 105, 60 

O 

609 

60 

DO 100 K=1,N 

D 

610 

70 

IF (IPIVOT(K)-l) 80, 100, 740 

O 

611 

80 

IF(A3S(AMAX)-ABS(A(J,K)) ) 85, 100, 100 

D 

612 

85 

IROW= J 

O 

613 

90 

ICOLUM=K 

D 

614 

95 

A MAX= A (J , K) 

D 

615 

100 

CONTINUE 

O 

616 

105 

CONTINUE 

O 

617 

110 

IP I VOT (ICOLUM) =IPlVOT < ICO LUM) *1 

D 

618 



O 

619 


INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

O 

620 



D 

621 

130 

IF (I ROW -ICO LUM) 140, 2*0, 140 

D 

622 

140 

CONTINUE 

O 

623 

150 

OO 200 L=1,N 

O 

624 

160 

SWAP=A(IROW,L> 

O 

625 

170 

A tIROW,L) =AC ICOLUMjL) 

D 

626 

200 

AdCOLUH ,L»=SWAP 

D 

627 

260 

INDEX(I,1) =IROW 

O 

628 

270 

INDEX (I,;2)=ICOLUM 

O 

629 

310 

PIVOT (I) = A (ICOLUM ,ISOLU1) 

0 

630 



0 

631 


DIVIOE PIVOT ROW BY PIVOT ELEMENT 

D 

632 



D 

633 

330 

A (ICOLUM , ICOLUM) =1.0 

D 

634 

340 

OO 350 L=1,W 

D 

635 

350 

A ( ICOLUM , L) = A (ICOLUM, L) ^P I VOT (I) 

O 

636 



D 

637 


REDUCE NON-P.IVOT ROWS 

O 

638 



D 

639 

380 

OO 550 L 1=1, N 

O 

640 

390 

IF (Ll-ICOLUM ) 400, 550, 400 

O 

641 

400 

T=A(L1, ICOLUM) 

O 

642 


109 



no o 


420 ACLi,ICOLUM) =0.0 
430 00 450 L=i*N 

450 A(L1,LI=A (L1,L>-A (ICOLUM, L)*T 
550 CONTINUE 

INTERCHANGE COLUMNS 

600 DO 710 1 = 1, N 
610 l = N«-i-I 

620 IF (IN0EX(L,1)-IN0EX(L,2) ) 630, 7 10, 
630 JROW=IND£X(L ,1) 

640 JC0LUM=INDEXi(L,:2) 

650 DO 705 K=1,N 
660 SWAP=A(K, JR3W) 

670 A <K , JROW ) = A( K ,'JC0 LUM) 

700 A(K, JCOLUM»=SHAP 
705 CONTINUE 
710 CONTINUE 
740 RETURN 
ENO 


D 643 
O 644 
O 645 
D 646 
D 647 
Q 648 
D 649 
D 650 
D 651 

630 D 652 

D 653 
O 654 
O 655 
O 656 
D 657 
D 658 
O 659 
O 660 
D 661 
D 662 
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SUBROUTINE RHS (X CPT , YC S T , ZC°T ,XW ,YW, ZW,DSM, GAMA M , SPAN, SPEEO, 
1GAMAK,NN,LL,NW,SINPHI,C0SPHI) 

C 

C THIS IS A SUBROUTINE TO COMPUTE THE RIGHT HAND SI OE OF THE MATRIX 
C EQUATION OEFINING THE VORTEX STRENGTHS. 

C 

OIMENSION Xlf(40), YW(40 ) , Z W (40 ) ,RH (2,2) , OS M(39 ) , V B AR (2 ) 

DIMENSION GA'MAK(10Q»1) 

OIMENSION SINPHK 20) , COSP HI ( 2 C) , X CPTC 14) , YCPT (10 ) ,ZCPT(10) 

P = 6.28 3 1 85 3 
M = 0 
C 

C CYCLE THROUGH CONTROL POINTS. 

900 00 50 1=1, NN 

DO 49 J=i,LL 
M = M «• 1 
C 

VYM = 0.0 
VZM = 0. 0 
SI NJ = SINPHI(J) 

COSJ = COSPHI(J) 

XCI = XCPT(I) 

YCJ = Y3PTU) 

ZCJ = ZCPT(J) 

C 

C COMPUTE VELOCITY INOUCEO BY MODEL. 

00 46 K= 1 , NW 
JJ = K 
00 45 L= 1,2 

RW(L,1) = SqRT((XW(JJ)-XCI)**2MYWfJJ>-YCJ)**2MZW(JJ)-ZCJ)**2) 
RW(L,2) = SQRT((XW(JJ)-KCI)**2+(YW(JJ>-YCJ)**2M7W(JJ>fZCJ> **2> 

JJ = KH 

43 CONTINUE 

00 44 L= 1 , 2 

44 VBAR(L) =-GAMAM*( DSM(K)** 2-(RW(l, L)-RW(2, L) > **2) * (RW(1,L> *RW(2,L> ) 
1/(P*RW(1 ,L)*RW(2,L)*(4.0* (RH(1,L) **2)*(0SM(K> **2>-{RW (1, L)**2-RW ( 2 
2,L)**2+0SM(K)**2) **2> ) 

L = K*1 

IF (COSJ. EQ. 0.0 GO TO 41 

VYM = V3AR(1)*((ZW(K)~ZCJ)*(XH(L) -XWTK) )- (XW( K) - XCI) * ( ?W ( L) -ZW(K) ) 
1)-VBAR(2) *((-ZW(K)-ZCJ)*< XW(L)-XH (K)I-(XW (K)-XCI)*(ZW (K)-ZM (L) )) 

2«- VYM 

41 IF (SINJ.EQ.0.0) GO TO 46 

VZM = (V3AR(1)-V3 AR(2) ) *( (XW(K)-XCI) * (YH( L)-YH(K) )- (YH(K)-YCJ)* 

1 (XW(L)-XH <K> ) ) «■ VZM 

45 CONTINUE 
C 

C STORE NORMAL VELOCITY IS GAMAK ARRAY ELEMENT M. 

54 GAMAK (M, 1) s V7H*SINJ - VYM*COSJ 

49 CONTINUE 

50 CONTINUE 
C 

C THE RIGHT HAND SIOE HAS BEEN GENERATED, RETURN TO CALLING PROGRAM. 

C 

RETURN 


0 663 

0 664 

D 665 
0 666 
0 667 

0 668 
D 669 
0 670 

0 671 

0 672 

0 673 

D 674 
0 675 

D 676 
0 677 

D 678 
0 679 

0 680 
D 631 
D 682 
D 633 
D 634 
D 685 
D 636 
0 687 

0 688 
0 689 

0 690 

D 691 
D 692 
0 693 

0 694 

0 695 

D 696 
0 697 

0 698 

0 699 

0 700 

0 701 

D 702 
0 703 

D 704 
0 705 

D 706 
D 707 
D 708 
0 709 

D 710 
0 711 

0 712 

D 713 
0 714 

D 715 
D 716 
0 717 


END 


D 718 


111 



SUBROUTINE WKIT ( XW, YW ,ZW ,X ,Y ,Z, S-INPHI ,COS*»HI ,SI OE, S, GAMA ,OSM, 0 

IGA MAN, SPEEO, SPAN, NW.NN.MM ,N1,LL, V Hi ,RHO,Q ,F AL ,F AD ,CHO RD, LIFT, DRAG , 0 
1STWK, VXTC.VYTC, ALPHAO, ALPHAI, ALFA A.VXMC.VYMC) 0 

C 0 

C THIS IS A SUBROUTINE TO ITERATE THE TRAILING VORTEX PAIR POSITION 0 

C ANO TO COMPUTE LIFT AND INDUCtO DRAG VALUES BASED UPON THE VELOCITY D 

C AT THE CENTER OF THE 80UN0 VORTEX. D 

C D 

DIMENSION X(15),Y (20) , Z< 2 0) , SINPH I ( 20) ,CO SPHI (20) .SIDE (20 ) ,S ( 14) 0 

DIMENSION GAMAI14.10) D 

DIMENSION XW (40) , YW(40),ZW(4O)»RW (Z,2)»DSM(39),VBAR(2) 0 

INTEGER A,B,C,0,E D 

LOGICAL STWK 0 

REAL LIFT D 

ALPHAI =0.0 D 

ALFAA s 0.0 0 

ALPHAO =0.0 D 

C 0 

C IF THIS IS TO BE THE LINEARIZED CASE. DO NOT ITERATE THE TRAILING 0 

C PAIR. GO DIRECTLY TO COMPUTE THE LIFT AND DRAG. 0 

IF (STHK) GO TO 704 0 

MMMM = NW-i 0 

C D 

C CYCLE THROUGH VERTICAL ANO LATERAL SHIFT OPERATIONS. 0 

00 701 LSHFT =1,2 D 

C 0 

C CYCLE THROUGH HAKE SEGMENTS. D 

DO 70 0 M = 2, MMMM D 

IF ((M.EQ.2) .ANO. (LSHFT. £0,2) ) GO TO 700 D 

C 0 

C SELECT COORDINATES FOR VEL02I TY COMPUTATION. NOTE ZCJ = 0.0 FOR CASE 0 
C OF FIRST TRAILING VORTEX SEGMENT. D 

XCI = XW(M) 0 

YCJ = YW (M) D 

IF (M.EQ.2) GO TO 20 D 

ZCJ = ZHCM) 0 

GO TO 33 D 

23 ZCJ =0.0 D 

30 CONTINUE D 

C D 

C COMPUTE VELOCITY AT THIS POINT. D 

CALL XYZVEL (XCI ,YC J,ZC J.XH ,YW, ZW ,X ,Y ,Z ,SI NPHI, COSPHI .SIDE, S 0 

1 .GAMA ,OSM,GAMAM,SPEEO,S»AN,NW,NN,MM,Ni,LL ,VXT,VYT,VZT .VXR.VYR, D 

1VZR.VXM, VYM, VZM) D 

VEL = SQRT (V XT**2 ♦ VYT**2 * V7T**2> D 

J = MU 0 

IF (M.NE.2) GO TO 743 D 

C D 

C IF THIS IS THE FIRST SEGMENT, COMPUTE NEW ANGLE OF ATTACK. 0 

ALPHAO = ASIN(”GAMAM*2.Z ( 6.283 IS 5 3* CHORD* VEL) ) 0 

ALPHAI = ATAN (VYT/VXT) 0 

ALFAA = ALPHAO ♦ ALPHAI D 

C 0 

C COMPUTE COORDINATE SHIFT. 0 

XSHFT = OSM(1)*COS(ALFAA) ♦ XW(1> - XW(2) 0 

YSHFT = QSM(i)*SIN( ALFAA) ♦ YW(1) - YW(2) 0 


719 

720 
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722 

723 
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725 

726 

727 

728 
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730 
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742 

743 
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746 

747 

748 

749 

750 

751 

752 
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755 
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759 

760 
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763 

764 

765 

766 

767 

768 

769 

770 

771 

772 

773 

774 
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ZSHFT = 0*0 
GO TO 57 
74 3 CONTINUE 

OCWX = VX T/VEL 
OCWY = VYT/VEL 
OCWZ = V7T/VEL 
XSHFT = DSM(M)*DCWX ♦ XH(M) 

XSHFT = XSHFIT - XW(J> 

IF (LSHFT.EQ.2) GO TO 149 
YSHFT = OSM(M)*DCWY «■ YU ( M) 

YSHFT = YSHFT - YW(J) 

GO TO 57 

149 ZSHFT = OSMC M) *DC WZ ♦ ZW(M) 

ZSHFT = ZSHFT - ZW(J) 

C 

C SHIFT ALL COORDINATES OOWNSTREAM OF VELOCITY COMPUTATION POINT. 

57 00 748 L=J,NW1 

K=L-1 

XH(L) = XH (L ) ♦ XSHFT 
IF (LSHFT.E9.2) GO TO 59 
59 YW(L) = YWCL) ♦ YSHFT 

GO TO 14 8 

59 ZW(L) = Z M (L ) ♦ ZSHFT 

C 

C COMPUTE NEW SEGMENT LENGTH. 

148 OSM(K) =' SQRTKXW <L> -X W( K ) ) ** 2«- ( Y W <L ) -YW( K) ) **2* (ZW (L ) -ZW (K) ) **2 ) 
74 8 CONTINUE 

790 CONTINUE 

701 CONTINUE 

4150 FORMAT ( 18H0WAKE COORDINATES ,/, 9X,’2HXW , 13X,2HVW,13X,2HZW) 

C 

C WRITE RESULT OF ITERATION PROCESS. 

WRITE (6.4150) 

4160 FORMAT (3F15.5) 

XCI = XW(1) 

DO 703 I = 1 , NW 

YCJ = XW(I) - XCI 

WRITE (6.4160) YC J, YW ( I) , ZW ( I ) 

70 3 CONTINUE 

704 CONTINUE 

C 

C COMPUTE LIFT AND INOUCEO DRAG OF WING, COMPARE WITH FREE AIR RESULT. 
XCI = XW (1) 

YCJ = YW(1) 

ZCJ =0. 

CALL XYVEL (XCI , YC J, ZC J.XW ,YW, ZW,X , Y.Z ,SI NPHI , COSPHI ,S IDE , S 

l.GAMA.OSM, GA'MAM, S PEED, SPAN, NW ,NN, NM.Nl.LL , VXT , VYT , VZT , VXR, VYR , 
1VZP,VXM, VYM, YZM) 

VXTC= VXT 
YYTC= VYT 
VXMC=VXM 
VYMC=VYN 

702 LIFT = RHO*SPAN*GAMAM 
ORAG * -LIFT*VYT 
LIFT = V XT*L I FT 

CLAR = ( (3.14159/4.)**2)/ (0*(SPAM**2) ) 


0 

0 

0 

0 

0 

D 

0 

0 

0 

D 

D 

D 

D 

0 

0 
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D 

D 

D 

0 

0 

D 

D 

0 

D 

D 

0 

D 

D 

0 

D 

D 

0 

0 

0 

0 

0 

0 

D 

D 

0 

D 

0 

D 

0 

0 

D 

D 

D 

0 

0 

0 

0 

0 

0 
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817 
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COIAR = ORAG*CLAR 
CLAR = LIFT'CLAR 
DELTAL = LIFT - FAL 
OELTAO = ORAG-FAO 
IF (STWK) ALFAA=Q • 0 
ALFAA=-A LFAA 
ALPHA0=“ ALPHA 0 

C 

C WRITE RESULTS OF COMPUTATIONS • 

WRITE (6,4176) ALFAA, ALPHAO, AL°HAI 

4176 FORMAT (2X,7HALFAA =, F7.)4 , 3X , 8HAL PH AO =,F7.4,'3X, 8HALPHAI =,F7.4) 
WRITE (5,4175) LIFT, DELTAL 

4175 FORMAT (12H3WING LIFT =,F 10. 5 ,4X,22HCHANGE DUE TO TUNNEL = ,F10.5> 
WRITE (6,4177) DRAG, DELTA D 

4177 FORMAT (12H0WING DRAG =, F 10. 5 ,4X, 22HCHANGE OUE TO TUNNEL =,F10.5> 
WRITE (6,4200) VXT,VYT 

4200 FORMAT ( 1H0, 39HT0TAL VELOCITIES AT WING CENTER VX = ,F1G.4,5X, 
15HVY = , F 10* 4 ) 

WRITE (6,4210) CLAR.CDIAR 

4210 FORMAT ( 1H0, 29HME ASUREO IN TUNNEL CL/AR = , F10 . 5, 5X ,9HCDI/AR = , 
1F10.5) 

RETURN TO CALLING PROGRAM. 

RETURN 
END 


0 

831 

0 

832 

D 

8 33 

0 

8 34 

0 

8 35 

0 

8 36 

0 

837 

0 

838 

D 

839 

D 
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9 
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0 
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D 
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D 
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0 
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0 
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0 
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0 
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0 
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D 
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D 
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0 
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D 
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D 
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SUBROUTINE XYZVEL (XCI ,Y C J ,ZC J,XW , YW. ZW ,X , Y , Z ,SI NPHI , COSPHI .SIDE ,S 0 857 


l,GAMA,OSM t GAMAM,SP££0,S»AN,NW,NN, MM.Nl.LL ,VXT,VYT,VZT,VXR,VYR, 0 858 

1VZR.VXM, VYM, VZM) 0 859 

C 0 860 

C THIS IS A SU 3 ROUTINE TO COMPUTE VELOCITY COMPONENTS INDUCED BY 0 861 

C TUNNEL AND LIFTING SYSTEM. 0 862 

C D 863 

DIMENSION XC15) ,Y (20) ,Z(20> ,SINPMI(20) ,CDSPHI(20) ,SI0E(20> ,St 14) , 0 864 

1R(15, 2G) , HL( 15,20 ) ,H0 (20 ) ,HYZ (20) ,XWC4Q> * YW (4 0 > , Z W (40 ) , RW (2 ,2) , D 865 

10SM(39>, VBAR(2),GAMA(14,10) D 866 

INTEGER A » 8, C » D.E D 867 

LOGICAL XONLY.XNY ,YNZ O 868 

C D 869 

C SET LOGICAL VARIABLES TO COMPUTE ONLY VELOCITY COMPONENTS REQUIRED. D 870 

IF (ZCJ. EQ.0.) GO TO 10 O 871 

XONLY = .FALSE. D 872 

XNY s .FALSE. D 873 

YNZ = .FALSE. O 874 

GO TO 643 O 875 

ENTRY XV EL D 876 

XONLY = .TRUE. D 8 77 

XNY = .FALSE. O 878 

YNZ = .FALSE. O 879 

GO TO 64 3 O 830 

ENTRY XYVEL D 881 

10 CONTINUE O 882 

XNY = .TRUE. D 883 

XONLY = .FALSE. O 834 

YNZ = .FALSE. O 885 

GO TO 64 3 O 8 36 

ENTRY YZVEL D 837 

YNZ = .TRUE. D 888 

XONLY = .FALSE. D 889 

XNY = .FALSE. D 890 

64 3 XTP = XCI O 891 

YTP = YCJ O 892 

ZTP = ZCJ O 893 

C D 894 

C COMPUTE LOCATION OF VORTEX MET WITH RESPECT TO POINT OF VELOCITY D 895 

C COMPUTATION. D 896 

DO 127 J s 1 , MM O 897 

HO ( J) = SQRC( (YTP-YCJ) >**2 + (ZTP - Z(J))**2) D 898 

HYZ(J) = SORT (( (ZTP-Z (J) ) *SINPHI { J>-{ YTP- Y(J) >*COSPHI (J) >**2) D 899 

DO 127 I = 1,N1 O 900 

R(I,J» = SORT ( (XTP-X(I)) **2 ♦ (YTP-YC J))**2 ♦ (ZTP-Z( Jl ) **21 D 901 

127 HLU,J>=SORT<(X(I)-XTP)**2 ♦ HYZ(J)**2) D 932 

VXR = 0. 60000 D 903 

VYR = 0. DC. COO O 904 

VZR = 0. 00000 0 905 

C O 906 

C CYCLE THROUGH VORTEX RECTANGLES. O 907 

DO 150 K - i,NN O 908 

00 150 L = 1 i LL 0 909 

C 0 910 

C SELECT VARIABLES' FOR THIS PARTICULAR VORTEX RECTANGLE OR RECTANGLES. 0 911 

A = L - 1 0 912 
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B = L 0 

C = LL*2 - L 0 

0 = C - 1 0 

E = K «■ 1 D 

IF (L - 1) 150, 129, 125 0 

125 IF (LL-L) 150,129,120 0 

128 RKA = R( K , A) 0 

RKC = R( K , C> > 0 

REA = R(E , A) 0 

REC = RIE,C> D 

HLKC = HL (K, C) 0 

HLEC = HL (E, C) 0 

HOA = HO (A) 0 

HOC = HD(C) 0 

YA = YCA) 0 

ZA = Z(A) 0 

ZC = ZO) 0 

HYZA = HYZ(A) 0 

HYZC = HYZ(C) 0 

129 SI NL = SINPHI(L) D 

COSL = COSPHI(L) 0 

RK8 = R(K,B) 0 

RKD = R(K,0) 0 

REfi = RfE.B) 0 

REO = R( E ,0) D 

HL KB = HL (K, B) 0 

HLEB = HL (E, 9) 0 

HOB = HO (B) 0 

HOO = H0(D) D 

SIDES * SIDE (B) 0 

OK = SCKI 0 

YB = Y(B) 0 

ZB = ZCS) 0 

ZO = Z(O) D 

XK = X(K) 0 

XE = X (E ) 0 

HYZB = HYZ(B) 0 

HYZO = HYZ(O) 0 

P = 25.13274 0 

0 

DETERMINE WHETHE R OR NOT YORTEX RECTANGLE LIES ON PLANE OF SYMMETRY. 0 
IF (L-l) 150, 130 , 131 0 

131 IF CLL-C) 132,130,150 0 

130 CONTINUE 0 

C 0 

C VORTEX RECTANGLE LIES ON PLANE OF SYMMETRY, USE FOLLOWING EQUATION TO 0 
C COMPUTE VELOCITY COMPONENTS TAKING SPECIAL CASES INTO ACCOUNT. 0 

VXPS =0.0 D 

VYPS =0.0 0 

VZPS = 0.0 D 

IF ( Y NT) GO TO 135 0 

VXPS = l./(P*SIDE BJMHYZB* ((RKDfRKB)* CSIDFB** 2- ( RKO-RKB)* *2)/ ( ( 0 

1HLKB**2> *PK0*RK9) - ( RE0+RE3) * (SI OEB**2-( RE0-RE8 > **2) / ( (HLEB**2) 0 

1* REO* RES ) ) )* GAMA( K,L) 0 

IF (XONLY) GO TO 72 0 

135 IF (COS. • EQ. 0.0) GO TO b 6 0 
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VYPS = (C0SL/(P*SIDE3)*(- ( (RKO*RK B) * ( SIDE B**2-(RKD-RKB) **2) /( ( 
2HLKB**2) *RK9*RK3) )*{XK-XTP> * ( (REO + REB) * (SI0E3* * 2 -(RED-REB)** 

22) / ( ( HLE 8**2 ) *REO *RE3) )* ( XE-XTP) ) * 1./ (P'*OK) * ( ( (RK3*RE8) *<QK** 

2 2- (RKB-RE 8)**2)/( (H03**2) *RK9*RE3 ) ) * ( ZB-Z TP) - ( (RKD*R£ O) * < OK** 2- 
2 (RKD-REO ) **2)/ ( <H 00**2 )*RKO* RED) ) *{ ZO-ZTP) )) *GAHA (K,L ) 

IF CX NY) GO TO 72 
GO TO 67 

65 VYPS = <1./(P*DK> •*(( (RK3*REB)*(0K** 

2 2- (RKB-RE 8)** 2) /( <H03**2> *PK 8*RE3 ) ) * ( ZB-ZTP) - ( ( RKO+REO) ♦ ( 0K**2- 
2 (RKO-REO) **2)/ ( (H00**2) *RKD*REO) ) *{ ZO-ZTP ))) *GAM A CK,L ) 

IF (XNY) GO TO 72 

67 VZPS = (i»/(P*OK)*(( (RK9* REO) *(0K ** 2 - ( RK D-RE D) **2) t { (H00**2> *RKO* 
3RE0) -RKB * RE 3) * (3K**2- (RKB-RE3 )**2)/{{HQ9**2) *RKB*REB) )* 
3CYB-YTP) ) )*GAMA(K,L) 

VZR = VZR * VZPS 
72 VXR = VXR * VXPS 
VYR = VYR * VYPS 
GO TO 150 
C 

C VORTEX RECTANGLES DO NOT LIE ON PLANE OF SYMMETRY, USE FOLLOWING 
C EQUATIONS TO COMPUTE VELOCITY COMPONENTS TAKING VARIOUS SPECIAL CASES 
C INTO ACCOUNT. 

13 2 CONTINUE 
VX = 0.0 
VY = 0.0 
VZ =0.0 

IF (Y NZ) GO TO 140 

VX =(i.Z (P*SI0E9> *((HYZ1* < (RKA+RKB) *(SI0E‘B**2 - ( PKA-RK9) **2) / C ( 
1HLK3**2> *RKA*RK3) - ( REA* REB) * CSI 0EB**2 - (REA-RE B> **2)/ (( HLE3**2) 
1*REA* REB ) ) ) * (HYZC*C <RKO*RKC)*(S I0E8**2 -(RKD-RKC) *♦ 2) / ( C 
1HLKC**2) *RKC*RKO) - (RE3+ REC) ♦ (SI 0EB**2 - (REO-REC) **2)/( (HLEC**2> 
1*REC*REB ) ))) ) *GAM A(K,L> 

IF (XONLY) GO TO 73 
140 IF (COSL.EO. 0.0) GO TO 58 

VY = (COSL/ (P*SIDE3)*(-((RKA*RK5)*( SIDE 3**2- (RKA-RKB) ** 2)/ ( ( 

2HL KB**2) *RK4*RK3) ♦ (RKO*RKC) *(SI DEB* *2 - (RKC-RKD)**2)/( (HLKC**2) 
2*RKC*RK0) )*( XK-XTP) ♦ ( ( REA*RE9) * (SIOEB** 2 - (REA-REB) **2) / ( ( 
2HLEB**2) * REA* REB) ♦ (REC* RED) *(SI OEB**2 - (REC -RED) **2)/((HLEC**2) 
2*REC*REO) )*(XE-XTP)) * 1./ (P*0K ) ♦ ( ( (RKB *RE3) * ( DK**2 - ( 

2RKB-REB) **2) / ( (HD 3**2 ) *RK B*RE B) ) * (ZB-ZTP) -( (RKO*RED)* (DK**2-(RKD- 
2RE0)**2) / ( (H00**2 )*RKD*RE 0) )* (ZD- ZTP) ♦ ( (R KC*REC) * (0K**2- (RKC-REC) 
2**2)/ ((HOC** 2) *RKC*REC)) * (ZC-ZTP) -( (RKA*REA) *:(0K**2- (RKA-RE A) **2) 
2/ ( (HO A** 2) *RKA*RE A) ) * ( ZA- ZTP) ) )*G AMA(K,L) 

IF (XNY) GO TO 73 
GO TO 69 

63 VY a ■ (i./(P*DK)*(( (RKB*REB)*(OK**2 -{ 


0 969 

0 970 

0 971 

D 972 
0 973 

0 974 

0 975 

D 976 
D 977 
D 978 
D 979 
0 930 

0 931 

0 982 

0 933 

D 984 
D 935 
D 986 
0 987 

0 988 

0 939 

0 990 

0 991 

D 992 
D 993 
0 994 

D 995 
0 996 

0 997 

0 998 

0 999 

D 1030 
D 1001 
D 1002 
0 1003 
D 1004 
0 1005 
0 1006 
0 1C 07 
0 1008 
D 1009 
0 1010 
D 1011 
0 1012 
0 1013 
0 1014 


2RKB-REB) * *2) / ( (HO 8** 2 ) *RK B*RE B) ) * (ZB-ZTP) -( (RKO* RED)* (OK**2-( RKD- 0 1015 
2REO)**2)/( (H 00**2 ) *RKO*RE D) ) * (ZD- ZTP) *((RKC*REC) * (DK**2- ( RKC- REC) D 1016 
2** 2)/ ( (HOC** 2) *RK C*REC) ) * ( ZC-ZTP) - ((RKA*REA)*(D<**2-( RKA-RE A) **2) D 1017 
2/ ( (HO A** 2 ) *RKA*RE A) ) * ( ZA- ZTP) ) ) *G AMA CK,L) D 1018 

IF (XNY) GO TO 73 0 1019 

69 IF (SINL.EQ. 0.00000) GO TO 70 0 1020 

VZ = (SI NL/( P*SIO EB) * ( ( ( RKA*RKB) * (SI0E9**2 -( RKA-RKB) **2) /( ( 0 1021 

3HLKB**2)*RKA*RK3) - (RKC*RKO) *(SI OEB**2 - (RKC-RKD) **2 ) / ( ( HL KC**2 ) 0 1022 

3*RKC*RKD> )*CXK-XTP) * ( ( REC+REO) * (SIOEB** 2 -(REC-RED) **2)/ ( ( 0 1023 

3HLEC**2) *REC*RE0) - (REA* REB) *(SI 0EB**2 - (REA-RE B)**2 )/(( HLE9**2) 0 1024 
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73 


71 

73 


3*REA*RE3) >*( XE-XTP)) * l./(P*OlO*(<{RKA+REA>*(DK**2 -crka 
3-REA)**2>/(<HDA**2>*RKA*R£A> - (R KOREC) * ( DK**2 - (RKC-REC)**2)/ ( ( 
3HDC**2)*RKC*REC>) *(YA-YFP ) ♦ ( C RKO* RED) * (0K**2 -CRKD-RED) ** 2) / 

3<(HDD**2>*RKD*RED> - (RK3 + RE9 )* CO K**2 - (RKB-REB) **2)/ ( <HDB**2> * 
3RK8*REB> ) *(YB-YTP ))) *GAMA <K,L) 

GO TO 71 

VZ = Cl./(P*OtO* < < (RKA+REA)* (DK**2 - (RKA 

3-REA)**2)/((!HDA**2)*RKA*REA) - (RKC+REC)* (0K**2 - CRKC-REC > **21/ < ( 
3HQC**2)*RKC*REC>) *(YA-Yr°) C { RKD* REO) * (DK**2 -<RKD-RE0) **2)/ 

3 < (HD0**2> *RK0*RED> - ( RK8 *R£B) *< Q K**2 - {R KB-RE3) **2) / C (HDB**2) * 
3RK8*R£B) ) * (Y 3-YTP )> > *G AMA <K,L) 

VZR = VZ R + VZ 


VXR a 
VYR = 


VXR 

VYR 


VX 

VY 


C 

C NOW 
15 0 


INDUCED BY MODEL, 


75 0 


■X5I )**2*(YM CJ)-YCJ) **2* tZW { J)-ZCJ) **2) 
-XCI ) **2* (YW (J)-YCJ) **2*<ZM ( J) *ZCJ)**2> 


74 5 


73 0 
744 


746 

C 


COMPUTE VELOCITY 
CONTINUE 
P = 6.2831853 
VXM = 0. 0 
VYM a 0. 0 
VZM a 0. 0 
DO 746 K=1,NW 
J = K 

DO 745 L=l,2 
RM(L,1) = SQRTUXW(J) 

RM(L,2) = SORT < (XW(J) 

J * K ♦ 1 
CONTINUE 
DO 744 L=l, 2 

H a 4.* <RM(1,LI**2> *<OSM (Kl **2> -!<RM(1,L> **2-RWt 2 ,L> **2*DSM <K> **2) 
1**2 

IF CH.LT. Ml.E-6) *4. *0SM < K) ** 2) ) GO TO 730 

VBAR(L) =-GAMAM*(0SN(K)**2-<RW<l, U-RM(2,L> >**2> * <RM(1,L> *RM(2,L> » 
1/(P*RW(1 ,L)*RW(2, L»*H) 

GO TO 744 
V3AR (L) = 0.0 
CONTINUE 
L a K*1 

IF <YNZ) GO TO 75 0 

VXM a VBARtil* l (Y M (K> - Y3 J > * < ZM<L) -ZWI K>> - CZM( K) -ZCJ) * (YM(L) -YM (K> ) 
1) - VfiAR (2)* ( (YM(K)-YCJ) * tZM < K)-Z MIL) > -<- ZMCKJ-ZC J) * CYM (L > -YM (K) > > 
2 * VXM 

IF CXONLY) GO TO 746 
CONTINUE 

VYM a V3AR(1) * 1 <Z M (K> -ZS J ) * ( XM (L) -XWCK) )- (XM(K)-XCI)* (ZW(L)-ZN(K) ) 
1)-VBAR(2> *(<-ZM<K)-ZCJ>MXM{L>-XM <K) ) - {XM (K) -XCI > * C ZM (K) -ZM (L ) ) ) 

2* VYM 

IF (XNY) GO TO 746 

VZM = IVBAR(1)-V3AR(2))*( (XM { K) -X Cl) * (YMC L) -YH(K ) )- (Y M f K) -YCJ) * 
lCXHCL)-XM(K) ) ) ♦ VZM 
CONTINUE 

VXT * VXM*VXR*SPEED 
VYT a VYM* VYR 
VZT a VZM*VZR 
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